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  // cluster removal loop
231  const T *firstOffset = &oldClusters.data().front();
232  for (typename DSV::const_iterator itdet = oldClusters.begin(), enddet = oldClusters.end(); itdet != enddet; ++itdet) {
233  DS oldDS = *itdet;
234 
235  if (oldDS.empty())
236  continue; // skip empty detsets
237 
238  uint32_t id = oldDS.detId();
239  DSF outds(*output, id);
240 
241  for (typename DS::const_iterator it = oldDS.begin(), ed = oldDS.end(); it != ed; ++it) {
242  uint32_t index = ((&*it) - firstOffset);
243  if (isGood[index]) {
244  outds.push_back(*it);
245  refs.push_back(index);
246  //std::cout << "HLTTrackClusterRemoverNew::cleanup " << typeid(T).name() << " reference " << index << " to " << (refs.size() - 1) << std::endl;
247  }
248  }
249  if (outds.empty())
250  outds.abort(); // not write in an empty DSV
251  }
252  // std::cout<<"fraction: "<<fraction<<std::endl;
253  if (oldRefs != nullptr)
254  mergeOld(refs, *oldRefs);
255  return output;
256 }
257 
258 void HLTTrackClusterRemoverNew::process(OmniClusterRef const &clusterReg, uint32_t subdet) {
259  if (clusterReg.id() != stripSourceProdID)
260  throw cms::Exception("Inconsistent Data")
261  << "HLTTrackClusterRemoverNew: strip cluster ref from Product ID = " << clusterReg.id()
262  << " does not match with source cluster collection (ID = " << stripSourceProdID << ")\n.";
263 
264  if (collectedRegStrips_.size() <= clusterReg.key()) {
265  edm::LogError("BadCollectionSize") << collectedRegStrips_.size() << " is smaller than " << clusterReg.key();
266 
267  assert(collectedRegStrips_.size() > clusterReg.key());
268  }
269  collectedRegStrips_[clusterReg.key()] = true;
270 }
271 
273  DetId detid = hit->geographicalId();
274  uint32_t subdet = detid.subdetId();
275 
276  assert((subdet > 0) && (subdet <= NumberOfParamBlocks));
277 
278  // chi2 cut
279  if (chi2 > pblocks_[subdet - 1].maxChi2_)
280  return;
281 
283  // std::cout<<"process pxl hit"<<std::endl;
284  if (!doPixel_)
285  return;
286  // this is a pixel, and i *know* it is
287  const SiPixelRecHit *pixelHit = static_cast<const SiPixelRecHit *>(hit);
288 
289  SiPixelRecHit::ClusterRef cluster = pixelHit->cluster();
290  if (cluster.id() != pixelSourceProdID)
291  throw cms::Exception("Inconsistent Data")
292  << "HLTTrackClusterRemoverNew: pixel cluster ref from Product ID = " << cluster.id()
293  << " does not match with source cluster collection (ID = " << pixelSourceProdID << ")\n.";
294 
295  assert(cluster.id() == pixelSourceProdID);
296  //DBG// cout << "HIT NEW PIXEL DETID = " << detid.rawId() << ", Cluster [ " << cluster.key().first << " / " << cluster.key().second << " ] " << endl;
297 
298  // if requested, cut on cluster size
299  if (pblocks_[subdet - 1].usesSize_ && (cluster->pixels().size() > pblocks_[subdet - 1].maxSize_))
300  return;
301 
302  // mark as used
303  //pixels[cluster.key()] = false;
304  assert(collectedPixels_.size() > cluster.key());
305  collectedPixels_[cluster.key()] = true;
306  } else { // aka Strip
307  if (!doStrip_)
308  return;
309  const type_info &hitType = typeid(*hit);
310  if (hitType == typeid(SiStripRecHit2D)) {
311  const SiStripRecHit2D *stripHit = static_cast<const SiStripRecHit2D *>(hit);
312  //DBG// cout << "Plain RecHit 2D: " << endl;
313  process(stripHit->omniClusterRef(), subdet);
314  // int clusCharge=0;
315  // for ( auto cAmp : stripHit->omniClusterRef().stripCluster().amplitudes() ) clusCharge+=cAmp;
316  // std::cout << "[HLTTrackClusterRemoverNew::process (SiStripRecHit2D) chi2: " << chi2 << " [" << subdet << " --> charge: " << clusCharge << "]" << std::endl;
317  } else if (hitType == typeid(SiStripRecHit1D)) {
318  const SiStripRecHit1D *hit1D = static_cast<const SiStripRecHit1D *>(hit);
319  process(hit1D->omniClusterRef(), subdet);
320  // int clusCharge=0;
321  // for ( auto cAmp : hit1D->omniClusterRef().stripCluster().amplitudes() ) clusCharge+=cAmp;
322  // std::cout << "[HLTTrackClusterRemoverNew::process (SiStripRecHit1D) chi2: " << chi2 << " [" << subdet << " --> charge: " << clusCharge << "]" << std::endl;
323  } else if (hitType == typeid(SiStripMatchedRecHit2D)) {
324  const SiStripMatchedRecHit2D *matchHit = static_cast<const SiStripMatchedRecHit2D *>(hit);
325  //DBG// cout << "Matched RecHit 2D: " << endl;
326  process(matchHit->monoClusterRef(), subdet);
327  // int clusCharge=0;
328  // for ( auto cAmp : matchHit->monoClusterRef().stripCluster().amplitudes() ) clusCharge+=cAmp;
329  // std::cout << "[HLTTrackClusterRemoverNew::process (SiStripMatchedRecHit2D:mono) chi2: " << chi2 << " [" << subdet << " --> charge: " << clusCharge << "]" << std::endl;
330 
331  process(matchHit->stereoClusterRef(), subdet);
332  // clusCharge=0;
333  // for ( auto cAmp : matchHit->stereoClusterRef().stripCluster().amplitudes() ) clusCharge+=cAmp;
334  // std::cout << "[HLTTrackClusterRemoverNew::process (SiStripMatchedRecHit2D:stereo) chi2: " << chi2 << " [" << subdet << " --> charge: " << clusCharge << "]" << std::endl;
335 
336  } else if (hitType == typeid(ProjectedSiStripRecHit2D)) {
337  const ProjectedSiStripRecHit2D *projHit = static_cast<const ProjectedSiStripRecHit2D *>(hit);
338  //DBG// cout << "Projected RecHit 2D: " << endl;
339  process(projHit->originalHit().omniClusterRef(), subdet);
340  // int clusCharge=0;
341  // for ( auto cAmp : projHit->originalHit().omniClusterRef().stripCluster().amplitudes() ) clusCharge+=cAmp;
342  // std::cout << "[HLTTrackClusterRemoverNew::process (ProjectedSiStripRecHit2D) chi2: " << chi2 << " [" << subdet << " --> charge: " << clusCharge << "]" << std::endl;
343  } else
344  throw cms::Exception("NOT IMPLEMENTED")
345  << "Don't know how to handle " << hitType.name() << " on detid " << detid.rawId() << "\n";
346  }
347 }
348 
349 /* Schematic picture of n-th step Iterative removal
350  * (that os removing clusters after n-th step tracking)
351  * clusters: [ C1 ] -> [ C2 ] -> ... -> [ Cn ] -> [ Cn + 1 ]
352  * ^ ^ ^--- OUTPUT "new" ID
353  * |-- before any removal |----- Source clusters
354  * |-- OUTPUT "old" ID |----- Hits in Traj. point here
355  * | \----- Old ClusterRemovalInfo "new" ID
356  * \-- Old ClusterRemovalInfo "old" ID
357  */
358 
360  ProductID pixelOldProdID, stripOldProdID;
361 
362  const auto &tgh = &iSetup.getData(tTrackerGeom_);
363 
365  if (doPixel_) {
366  iEvent.getByToken(pixelClusters_, pixelClusters);
368  }
369 
371  if (doStrip_) {
372  iEvent.getByToken(stripClusters_, stripClusters);
374  }
375 
376  //Handle<TrajTrackAssociationCollection> trajectories;
378  iEvent.getByToken(trajectories_, trajectories);
379 
380  if (mergeOld_) {
383  iEvent.getByToken(oldPxlMaskToken_, oldPxlMask);
384  iEvent.getByToken(oldStrMaskToken_, oldStrMask);
385  LogDebug("TrackClusterRemover") << "to merge in, " << oldStrMask->size() << " strp and " << oldPxlMask->size()
386  << " pxl";
387  oldStrMask->copyMaskTo(collectedRegStrips_);
388  oldPxlMask->copyMaskTo(collectedPixels_);
389  collectedRegStrips_.resize(stripClusters->dataSize(), false);
390  } else {
391  collectedRegStrips_.resize(stripClusters->dataSize(), false);
392  collectedPixels_.resize(pixelClusters->dataSize(), false);
393  }
394 
395  //for (TrajTrackAssociationCollection::const_iterator it = trajectories->begin(), ed = trajectories->end(); it != ed; ++it) {
396  // const Trajectory &tj = * it->key;
397 
398  for (std::vector<Trajectory>::const_iterator it = trajectories->begin(), ed = trajectories->end(); it != ed; ++it) {
399  const Trajectory &tj = *it;
400  const std::vector<TrajectoryMeasurement> &tms = tj.measurements();
401 
402  std::vector<TrajectoryMeasurement>::const_iterator itm, endtm;
403  for (itm = tms.begin(), endtm = tms.end(); itm != endtm; ++itm) {
404  const TrackingRecHit *hit = itm->recHit()->hit();
405  if (!hit->isValid())
406  continue;
407  // std::cout<<"process hit"<<std::endl;
408  process(hit, itm->estimate(), tgh);
409  }
410  }
411 
412  // std::cout << " => collectedRegStrips_: " << collectedRegStrips_.size() << std::endl;
413  // std::cout << " total strip to skip (before charge check): "<<std::count(collectedRegStrips_.begin(),collectedRegStrips_.end(),true) << std::endl;
414 
415  // checks only cluster already found! creates regression
416  if (doStripChargeCheck_) {
417  // std::cout << "[HLTTrackClusterRemoverNew::produce] doStripChargeCheck_: " << (doStripChargeCheck_ ? "true" : "false") << " stripClusters: " << stripClusters->size() << std::endl;
418 
419  auto const &clusters = stripClusters->data();
420  for (auto const &item : stripClusters->ids()) {
421  if (!item.isValid())
422  continue; // not umpacked
423 
424  DetId detid = item.id;
425  uint32_t subdet = detid.subdetId();
426  if (!pblocks_[subdet - 1].cutOnStripCharge_)
427  continue;
428 
429  // std::cout << " i: " << i << " --> detid: " << detid << " --> subdet: " << subdet << std::endl;
430 
431  for (int i = item.offset; i < item.offset + int(item.size); ++i) {
432  int clusCharge = 0;
433  for (auto cAmp : clusters[i].amplitudes())
434  clusCharge += cAmp;
435 
436  // if (clusCharge < pblocks_[subdet-1].minGoodStripCharge_) std::cout << " clusCharge: " << clusCharge << std::endl;
437  if (clusCharge < pblocks_[subdet - 1].minGoodStripCharge_)
438  collectedRegStrips_[i] = true; // (|= does not work!)
439  }
440  }
441  }
442 
443  // std::cout << " => collectedRegStrips_: " << collectedRegStrips_.size() << std::endl;
444  // std::cout << " total strip to skip: "<<std::count(collectedRegStrips_.begin(),collectedRegStrips_.end(),true) << std::endl;
445  LogDebug("TrackClusterRemover") << "total strip to skip: "
446  << std::count(collectedRegStrips_.begin(), collectedRegStrips_.end(), true);
447  // std::cout << "TrackClusterRemover" <<" total strip to skip: "<<std::count(collectedRegStrips_.begin(),collectedRegStrips_.end(),true)<<std::endl;
448  iEvent.put(std::make_unique<StripMaskContainer>(edm::RefProd<edmNew::DetSetVector<SiStripCluster> >(stripClusters),
450 
451  LogDebug("TrackClusterRemover") << "total pxl to skip: "
452  << std::count(collectedPixels_.begin(), collectedPixels_.end(), true);
453  iEvent.put(std::make_unique<PixelMaskContainer>(edm::RefProd<edmNew::DetSetVector<SiPixelCluster> >(pixelClusters),
455 
456  collectedRegStrips_.clear();
457  collectedPixels_.clear();
458 }
459 
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:307
ProductID id() const
Accessor for product ID.
Definition: Ref.h:238
T const & getData(const ESGetToken< T, R > &iToken) const noexcept(false)
Definition: EventSetup.h:119
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:244
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
std::vector< uint32_t > Indices
Definition: DetId.h:17
const_iterator begin(bool update=false) const
size_type size() const
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)
Definition: output.py:1
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)