CMS 3D CMS Logo

HLTTrackClusterRemoverNew.cc
Go to the documentation of this file.
6 
21 
24 
28 
33 
34 //
35 // class decleration
36 //
37 
39 public:
41  ~HLTTrackClusterRemoverNew() override;
42  void produce(edm::Event &iEvent, const edm::EventSetup &iSetup) override;
43 
44 private:
46  struct ParamBlock {
48  ParamBlock(const edm::ParameterSet &iConfig)
49  : isSet_(true),
50  usesCharge_(iConfig.exists("maxCharge")),
51  usesSize_(iConfig.exists("maxSize")),
52  cutOnPixelCharge_(iConfig.exists("minGoodPixelCharge")),
53  cutOnStripCharge_(iConfig.exists("minGoodStripCharge")),
54  maxChi2_(iConfig.getParameter<double>("maxChi2")),
55  maxCharge_(usesCharge_ ? iConfig.getParameter<double>("maxCharge") : 0),
56  minGoodPixelCharge_(cutOnPixelCharge_ ? iConfig.getParameter<double>("minGoodPixelCharge") : 0),
57  minGoodStripCharge_(cutOnStripCharge_ ? iConfig.getParameter<double>("minGoodStripCharge") : 0),
58  maxSize_(usesSize_ ? iConfig.getParameter<uint32_t>("maxSize") : 0) {}
61  size_t maxSize_;
62  };
63  static const unsigned int NumberOfParamBlocks = 6;
64 
65  bool doTracks_;
67  bool mergeOld_;
68 
76 
78  void readPSet(const edm::ParameterSet &iConfig,
79  const std::string &name,
80  int id1 = -1,
81  int id2 = -1,
82  int id3 = -1,
83  int id4 = -1,
84  int id5 = -1,
85  int id6 = -1);
86 
87  std::vector<uint8_t> pixels, strips; // avoid unneed alloc/dealloc of this
88  edm::ProductID pixelSourceProdID, stripSourceProdID; // ProdIDs refs must point to (for consistency tests)
89 
90  inline void process(const TrackingRecHit *hit, float chi2, const TrackerGeometry *tg);
91  inline void process(const OmniClusterRef &cluRef, uint32_t subdet);
92 
93  template <typename T>
94  std::unique_ptr<edmNew::DetSetVector<T> > cleanup(const edmNew::DetSetVector<T> &oldClusters,
95  const std::vector<uint8_t> &isGood,
97  const reco::ClusterRemovalInfo::Indices *oldRefs);
98 
99  // Carries in full removal info about a given det from oldRefs
101 
104  std::vector<bool> collectedRegStrips_;
105  std::vector<bool> collectedPixels_;
106 };
107 
108 using namespace std;
109 using namespace edm;
110 using namespace reco;
111 
113  const edm::ParameterSet &iConfig, const std::string &name, int id1, int id2, int id3, int id4, int id5, int id6) {
114  if (iConfig.exists(name)) {
115  ParamBlock pblock(iConfig.getParameter<ParameterSet>(name));
116  if (id1 == -1) {
117  fill(pblocks_, pblocks_ + NumberOfParamBlocks, pblock);
118  } else {
119  pblocks_[id1] = pblock;
120  if (id2 != -1)
121  pblocks_[id2] = pblock;
122  if (id3 != -1)
123  pblocks_[id3] = pblock;
124  if (id4 != -1)
125  pblocks_[id4] = pblock;
126  if (id5 != -1)
127  pblocks_[id5] = pblock;
128  if (id6 != -1)
129  pblocks_[id6] = pblock;
130  }
131  }
132 }
133 
136  doTracks_(iConfig.exists("trajectories")),
137  doStrip_(iConfig.existsAs<bool>("doStrip") ? iConfig.getParameter<bool>("doStrip") : true),
138  doPixel_(iConfig.existsAs<bool>("doPixel") ? iConfig.getParameter<bool>("doPixel") : true),
139  mergeOld_(false),
140  makeProducts_(true),
141  doStripChargeCheck_(
142  iConfig.existsAs<bool>("doStripChargeCheck") ? iConfig.getParameter<bool>("doStripChargeCheck") : false),
143  doPixelChargeCheck_(
144  iConfig.existsAs<bool>("doPixelChargeCheck") ? iConfig.getParameter<bool>("doPixelChargeCheck") : false)
145 
146 {
147  if (iConfig.exists("oldClusterRemovalInfo")) {
148  oldPxlMaskToken_ = consumes<PixelMaskContainer>(iConfig.getParameter<InputTag>("oldClusterRemovalInfo"));
149  oldStrMaskToken_ = consumes<StripMaskContainer>(iConfig.getParameter<InputTag>("oldClusterRemovalInfo"));
150  if (not(iConfig.getParameter<InputTag>("oldClusterRemovalInfo") == edm::InputTag()))
151  mergeOld_ = true;
152  }
153 
155  throw cms::Exception("Configuration Error")
156  << "HLTTrackClusterRemoverNew: Charge check asked without cluster collection ";
158  throw cms::Exception("Configuration Error")
159  << "HLTTrackClusterRemoverNew: Pixel cluster charge check not yet implemented";
160 
162  readPSet(iConfig, "Common", -1);
163  if (doPixel_) {
164  readPSet(iConfig, "Pixel", 0, 1);
165  readPSet(iConfig, "PXB", 0);
166  readPSet(iConfig, "PXE", 1);
167  }
168  if (doStrip_) {
169  readPSet(iConfig, "Strip", 2, 3, 4, 5);
170  readPSet(iConfig, "StripInner", 2, 3);
171  readPSet(iConfig, "StripOuter", 4, 5);
172  readPSet(iConfig, "TIB", 2);
173  readPSet(iConfig, "TID", 3);
174  readPSet(iConfig, "TOB", 4);
175  readPSet(iConfig, "TEC", 5);
176  }
177 
178  bool usingCharge = false;
179  for (size_t i = 0; i < NumberOfParamBlocks; ++i) {
180  if (!pblocks_[i].isSet_)
181  throw cms::Exception("Configuration Error")
182  << "HLTTrackClusterRemoverNew: Missing configuration for detector with subDetID = " << (i + 1);
183  if (pblocks_[i].usesCharge_ && !usingCharge) {
184  throw cms::Exception("Configuration Error")
185  << "HLTTrackClusterRemoverNew: Configuration for subDetID = " << (i + 1)
186  << " uses cluster charge, which is not enabled.";
187  }
188  }
189 
190  // trajectories_ = consumes<vector<Trajectory> >(iConfig.getParameter<InputTag>("trajectories"));
191  if (doTracks_)
192  trajectories_ = consumes<vector<Trajectory> >(iConfig.getParameter<InputTag>("trajectories"));
193  if (doPixel_)
194  pixelClusters_ = consumes<edmNew::DetSetVector<SiPixelCluster> >(iConfig.getParameter<InputTag>("pixelClusters"));
195  if (doStrip_)
196  stripClusters_ = consumes<edmNew::DetSetVector<SiStripCluster> >(iConfig.getParameter<InputTag>("stripClusters"));
197  if (mergeOld_) {
198  oldPxlMaskToken_ = consumes<PixelMaskContainer>(iConfig.getParameter<InputTag>("oldClusterRemovalInfo"));
199  oldStrMaskToken_ = consumes<StripMaskContainer>(iConfig.getParameter<InputTag>("oldClusterRemovalInfo"));
200  }
201 
202  //produces<edmNew::DetSetVector<SiPixelClusterRefNew> >();
203  //produces<edmNew::DetSetVector<SiStripRecHit1D::ClusterRegionalRef> >();
204 
205  produces<edm::ContainerMask<edmNew::DetSetVector<SiPixelCluster> > >();
206  produces<edm::ContainerMask<edmNew::DetSetVector<SiStripCluster> > >();
207 }
208 
210 
212  const ClusterRemovalInfo::Indices &oldRefs) {
213  for (size_t i = 0, n = refs.size(); i < n; ++i) {
214  refs[i] = oldRefs[refs[i]];
215  }
216 }
217 
218 template <typename T>
219 std::unique_ptr<edmNew::DetSetVector<T> > HLTTrackClusterRemoverNew::cleanup(
220  const edmNew::DetSetVector<T> &oldClusters,
221  const std::vector<uint8_t> &isGood,
223  const reco::ClusterRemovalInfo::Indices *oldRefs) {
224  typedef typename edmNew::DetSetVector<T> DSV;
225  typedef typename edmNew::DetSetVector<T>::FastFiller DSF;
226  typedef typename edmNew::DetSet<T> DS;
227  auto output = std::make_unique<DSV>();
228  output->reserve(oldClusters.size(), oldClusters.dataSize());
229 
230  unsigned int countOld = 0;
231  unsigned int countNew = 0;
232 
233  // cluster removal loop
234  const T *firstOffset = &oldClusters.data().front();
235  for (typename DSV::const_iterator itdet = oldClusters.begin(), enddet = oldClusters.end(); itdet != enddet; ++itdet) {
236  DS oldDS = *itdet;
237 
238  if (oldDS.empty())
239  continue; // skip empty detsets
240 
241  uint32_t id = oldDS.detId();
242  DSF outds(*output, id);
243 
244  for (typename DS::const_iterator it = oldDS.begin(), ed = oldDS.end(); it != ed; ++it) {
245  uint32_t index = ((&*it) - firstOffset);
246  countOld++;
247  if (isGood[index]) {
248  outds.push_back(*it);
249  countNew++;
250  refs.push_back(index);
251  //std::cout << "HLTTrackClusterRemoverNew::cleanup " << typeid(T).name() << " reference " << index << " to " << (refs.size() - 1) << std::endl;
252  }
253  }
254  if (outds.empty())
255  outds.abort(); // not write in an empty DSV
256  }
257  // double fraction = countNew / (double) countOld;
258  // std::cout<<"fraction: "<<fraction<<std::endl;
259  if (oldRefs != nullptr)
260  mergeOld(refs, *oldRefs);
261  return output;
262 }
263 
264 void HLTTrackClusterRemoverNew::process(OmniClusterRef const &clusterReg, uint32_t subdet) {
265  if (clusterReg.id() != stripSourceProdID)
266  throw cms::Exception("Inconsistent Data")
267  << "HLTTrackClusterRemoverNew: strip cluster ref from Product ID = " << clusterReg.id()
268  << " does not match with source cluster collection (ID = " << stripSourceProdID << ")\n.";
269 
270  if (collectedRegStrips_.size() <= clusterReg.key()) {
271  edm::LogError("BadCollectionSize") << collectedRegStrips_.size() << " is smaller than " << clusterReg.key();
272 
273  assert(collectedRegStrips_.size() > clusterReg.key());
274  }
275  collectedRegStrips_[clusterReg.key()] = true;
276 }
277 
279  DetId detid = hit->geographicalId();
280  uint32_t subdet = detid.subdetId();
281 
282  assert((subdet > 0) && (subdet <= NumberOfParamBlocks));
283 
284  // chi2 cut
285  if (chi2 > pblocks_[subdet - 1].maxChi2_)
286  return;
287 
289  // std::cout<<"process pxl hit"<<std::endl;
290  if (!doPixel_)
291  return;
292  // this is a pixel, and i *know* it is
293  const SiPixelRecHit *pixelHit = static_cast<const SiPixelRecHit *>(hit);
294 
295  SiPixelRecHit::ClusterRef cluster = pixelHit->cluster();
296  if (cluster.id() != pixelSourceProdID)
297  throw cms::Exception("Inconsistent Data")
298  << "HLTTrackClusterRemoverNew: pixel cluster ref from Product ID = " << cluster.id()
299  << " does not match with source cluster collection (ID = " << pixelSourceProdID << ")\n.";
300 
301  assert(cluster.id() == pixelSourceProdID);
302  //DBG// cout << "HIT NEW PIXEL DETID = " << detid.rawId() << ", Cluster [ " << cluster.key().first << " / " << cluster.key().second << " ] " << endl;
303 
304  // if requested, cut on cluster size
305  if (pblocks_[subdet - 1].usesSize_ && (cluster->pixels().size() > pblocks_[subdet - 1].maxSize_))
306  return;
307 
308  // mark as used
309  //pixels[cluster.key()] = false;
310  assert(collectedPixels_.size() > cluster.key());
311  collectedPixels_[cluster.key()] = true;
312  } else { // aka Strip
313  if (!doStrip_)
314  return;
315  const type_info &hitType = typeid(*hit);
316  if (hitType == typeid(SiStripRecHit2D)) {
317  const SiStripRecHit2D *stripHit = static_cast<const SiStripRecHit2D *>(hit);
318  //DBG// cout << "Plain RecHit 2D: " << endl;
319  process(stripHit->omniClusterRef(), subdet);
320  // int clusCharge=0;
321  // for ( auto cAmp : stripHit->omniClusterRef().stripCluster().amplitudes() ) clusCharge+=cAmp;
322  // std::cout << "[HLTTrackClusterRemoverNew::process (SiStripRecHit2D) chi2: " << chi2 << " [" << subdet << " --> charge: " << clusCharge << "]" << std::endl;
323  } else if (hitType == typeid(SiStripRecHit1D)) {
324  const SiStripRecHit1D *hit1D = static_cast<const SiStripRecHit1D *>(hit);
325  process(hit1D->omniClusterRef(), subdet);
326  // int clusCharge=0;
327  // for ( auto cAmp : hit1D->omniClusterRef().stripCluster().amplitudes() ) clusCharge+=cAmp;
328  // std::cout << "[HLTTrackClusterRemoverNew::process (SiStripRecHit1D) chi2: " << chi2 << " [" << subdet << " --> charge: " << clusCharge << "]" << std::endl;
329  } else if (hitType == typeid(SiStripMatchedRecHit2D)) {
330  const SiStripMatchedRecHit2D *matchHit = static_cast<const SiStripMatchedRecHit2D *>(hit);
331  //DBG// cout << "Matched RecHit 2D: " << endl;
332  process(matchHit->monoClusterRef(), subdet);
333  // int clusCharge=0;
334  // for ( auto cAmp : matchHit->monoClusterRef().stripCluster().amplitudes() ) clusCharge+=cAmp;
335  // std::cout << "[HLTTrackClusterRemoverNew::process (SiStripMatchedRecHit2D:mono) chi2: " << chi2 << " [" << subdet << " --> charge: " << clusCharge << "]" << std::endl;
336 
337  process(matchHit->stereoClusterRef(), subdet);
338  // clusCharge=0;
339  // for ( auto cAmp : matchHit->stereoClusterRef().stripCluster().amplitudes() ) clusCharge+=cAmp;
340  // std::cout << "[HLTTrackClusterRemoverNew::process (SiStripMatchedRecHit2D:stereo) chi2: " << chi2 << " [" << subdet << " --> charge: " << clusCharge << "]" << std::endl;
341 
342  } else if (hitType == typeid(ProjectedSiStripRecHit2D)) {
343  const ProjectedSiStripRecHit2D *projHit = static_cast<const ProjectedSiStripRecHit2D *>(hit);
344  //DBG// cout << "Projected RecHit 2D: " << endl;
345  process(projHit->originalHit().omniClusterRef(), subdet);
346  // int clusCharge=0;
347  // for ( auto cAmp : projHit->originalHit().omniClusterRef().stripCluster().amplitudes() ) clusCharge+=cAmp;
348  // std::cout << "[HLTTrackClusterRemoverNew::process (ProjectedSiStripRecHit2D) chi2: " << chi2 << " [" << subdet << " --> charge: " << clusCharge << "]" << std::endl;
349  } else
350  throw cms::Exception("NOT IMPLEMENTED")
351  << "Don't know how to handle " << hitType.name() << " on detid " << detid.rawId() << "\n";
352  }
353 }
354 
355 /* Schematic picture of n-th step Iterative removal
356  * (that os removing clusters after n-th step tracking)
357  * clusters: [ C1 ] -> [ C2 ] -> ... -> [ Cn ] -> [ Cn + 1 ]
358  * ^ ^ ^--- OUTPUT "new" ID
359  * |-- before any removal |----- Source clusters
360  * |-- OUTPUT "old" ID |----- Hits in Traj. point here
361  * | \----- Old ClusterRemovalInfo "new" ID
362  * \-- Old ClusterRemovalInfo "old" ID
363  */
364 
366  ProductID pixelOldProdID, stripOldProdID;
367 
368  const auto &tgh = &iSetup.getData(tTrackerGeom_);
369 
371  if (doPixel_) {
372  iEvent.getByToken(pixelClusters_, pixelClusters);
374  }
375 
377  if (doStrip_) {
378  iEvent.getByToken(stripClusters_, stripClusters);
380  }
381 
382  //Handle<TrajTrackAssociationCollection> trajectories;
384  iEvent.getByToken(trajectories_, trajectories);
385 
386  if (mergeOld_) {
389  iEvent.getByToken(oldPxlMaskToken_, oldPxlMask);
390  iEvent.getByToken(oldStrMaskToken_, oldStrMask);
391  LogDebug("TrackClusterRemover") << "to merge in, " << oldStrMask->size() << " strp and " << oldPxlMask->size()
392  << " pxl";
393  oldStrMask->copyMaskTo(collectedRegStrips_);
394  oldPxlMask->copyMaskTo(collectedPixels_);
395  collectedRegStrips_.resize(stripClusters->dataSize(), false);
396  } else {
397  collectedRegStrips_.resize(stripClusters->dataSize(), false);
398  collectedPixels_.resize(pixelClusters->dataSize(), false);
399  }
400 
401  //for (TrajTrackAssociationCollection::const_iterator it = trajectories->begin(), ed = trajectories->end(); it != ed; ++it) {
402  // const Trajectory &tj = * it->key;
403 
404  for (std::vector<Trajectory>::const_iterator it = trajectories->begin(), ed = trajectories->end(); it != ed; ++it) {
405  const Trajectory &tj = *it;
406  const std::vector<TrajectoryMeasurement> &tms = tj.measurements();
407 
408  std::vector<TrajectoryMeasurement>::const_iterator itm, endtm;
409  for (itm = tms.begin(), endtm = tms.end(); itm != endtm; ++itm) {
410  const TrackingRecHit *hit = itm->recHit()->hit();
411  if (!hit->isValid())
412  continue;
413  // std::cout<<"process hit"<<std::endl;
414  process(hit, itm->estimate(), tgh);
415  }
416  }
417 
418  // std::cout << " => collectedRegStrips_: " << collectedRegStrips_.size() << std::endl;
419  // std::cout << " total strip to skip (before charge check): "<<std::count(collectedRegStrips_.begin(),collectedRegStrips_.end(),true) << std::endl;
420 
421  // checks only cluster already found! creates regression
422  if (doStripChargeCheck_) {
423  // std::cout << "[HLTTrackClusterRemoverNew::produce] doStripChargeCheck_: " << (doStripChargeCheck_ ? "true" : "false") << " stripClusters: " << stripClusters->size() << std::endl;
424 
425  auto const &clusters = stripClusters->data();
426  for (auto const &item : stripClusters->ids()) {
427  if (!item.isValid())
428  continue; // not umpacked
429 
430  DetId detid = item.id;
431  uint32_t subdet = detid.subdetId();
432  if (!pblocks_[subdet - 1].cutOnStripCharge_)
433  continue;
434 
435  // std::cout << " i: " << i << " --> detid: " << detid << " --> subdet: " << subdet << std::endl;
436 
437  for (int i = item.offset; i < item.offset + int(item.size); ++i) {
438  int clusCharge = 0;
439  for (auto cAmp : clusters[i].amplitudes())
440  clusCharge += cAmp;
441 
442  // if (clusCharge < pblocks_[subdet-1].minGoodStripCharge_) std::cout << " clusCharge: " << clusCharge << std::endl;
443  if (clusCharge < pblocks_[subdet - 1].minGoodStripCharge_)
444  collectedRegStrips_[i] = true; // (|= does not work!)
445  }
446  }
447  }
448 
449  // std::cout << " => collectedRegStrips_: " << collectedRegStrips_.size() << std::endl;
450  // std::cout << " total strip to skip: "<<std::count(collectedRegStrips_.begin(),collectedRegStrips_.end(),true) << std::endl;
451  LogDebug("TrackClusterRemover") << "total strip to skip: "
452  << std::count(collectedRegStrips_.begin(), collectedRegStrips_.end(), true);
453  // std::cout << "TrackClusterRemover" <<" total strip to skip: "<<std::count(collectedRegStrips_.begin(),collectedRegStrips_.end(),true)<<std::endl;
454  iEvent.put(std::make_unique<StripMaskContainer>(edm::RefProd<edmNew::DetSetVector<SiStripCluster> >(stripClusters),
456 
457  LogDebug("TrackClusterRemover") << "total pxl to skip: "
458  << std::count(collectedPixels_.begin(), collectedPixels_.end(), true);
459  iEvent.put(std::make_unique<PixelMaskContainer>(edm::RefProd<edmNew::DetSetVector<SiPixelCluster> >(pixelClusters),
461 
462  collectedRegStrips_.clear();
463  collectedPixels_.clear();
464 }
465 
ClusterRef cluster() const
Definition: SiPixelRecHit.h:47
edm::EDGetTokenT< edmNew::DetSetVector< SiStripCluster > > stripClusters_
void copyMaskTo(std::vector< bool > &) const
Definition: ContainerMask.h:87
size_type dataSize() const
ESGetTokenH3DDVariant esConsumes(std::string const &Record, edm::ConsumesCollector &)
Definition: DeDxTools.cc:283
T getParameter(std::string const &) const
Definition: ParameterSet.h:303
ProductID id() const
Accessor for product ID.
Definition: Ref.h:244
static const unsigned int NumberOfParamBlocks
edm::EDGetTokenT< edmNew::DetSetVector< SiPixelCluster > > pixelClusters_
bool exists(std::string const &parameterName) const
checks if a parameter exists
OmniClusterRef const & stereoClusterRef() const
OmniClusterRef const & omniClusterRef() const
data_type const * data(size_t cell) const
Log< level::Error, false > LogError
assert(be >=bs)
OmniClusterRef const & monoClusterRef() const
key_type key() const
Accessor for product key.
Definition: Ref.h:250
edm::EDGetTokenT< PixelMaskContainer > oldPxlMaskToken_
edm::ProductID id() const
HLTTrackClusterRemoverNew(const edm::ParameterSet &iConfig)
const_iterator end(bool update=false) const
edm::ContainerMask< edmNew::DetSetVector< SiStripCluster > > StripMaskContainer
std::unique_ptr< edmNew::DetSetVector< T > > cleanup(const edmNew::DetSetVector< T > &oldClusters, const std::vector< uint8_t > &isGood, reco::ClusterRemovalInfo::Indices &refs, const reco::ClusterRemovalInfo::Indices *oldRefs)
DataContainer const & measurements() const
Definition: Trajectory.h:178
void readPSet(const edm::ParameterSet &iConfig, const std::string &name, int id1=-1, int id2=-1, int id3=-1, int id4=-1, int id5=-1, int id6=-1)
int iEvent
Definition: GenABIO.cc:224
void process(const TrackingRecHit *hit, float chi2, const TrackerGeometry *tg)
edm::ESGetToken< TrackerGeometry, TrackerDigiGeometryRecord > const tTrackerGeom_
ParamBlock pblocks_[NumberOfParamBlocks]
edm::EDGetTokenT< std::vector< Trajectory > > trajectories_
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:16
constexpr int subdetId() const
get the contents of the subdetector field (not cast into any detector&#39;s numbering enum) ...
Definition: DetId.h:48
bool getData(T &iHolder) const
Definition: EventSetup.h:122
std::vector< uint32_t > Indices
Definition: DetId.h:17
const_iterator begin(bool update=false) const
size_type size() const
constexpr uint32_t rawId() const
get the raw id
Definition: DetId.h:57
void mergeOld(reco::ClusterRemovalInfo::Indices &refs, const reco::ClusterRemovalInfo::Indices &oldRefs)
fixed size matrix
HLT enums.
const GeomDetEnumerators::SubDetector geomDetSubDetector(int subdet) const
SiStripRecHit2D originalHit() const
void produce(edm::Event &iEvent, const edm::EventSetup &iSetup) override
ParamBlock(const edm::ParameterSet &iConfig)
unsigned int key() const
size_t size() const
Definition: ContainerMask.h:53
edm::ContainerMask< edmNew::DetSetVector< SiPixelCluster > > PixelMaskContainer
long double T
bool isTrackerPixel(GeomDetEnumerators::SubDetector m)
edm::EDGetTokenT< StripMaskContainer > oldStrMaskToken_
Our base class.
Definition: SiPixelRecHit.h:23
#define LogDebug(id)