CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
TrackListMerger.cc
Go to the documentation of this file.
1 //
2 // Package: RecoTracker/FinalTrackSelectors
3 // Class: TrackListMerger
4 //
5 // Description: Hit Dumper
6 //
7 // Original Author: Steve Wagner, stevew@pizero.colorado.edu
8 // Created: Sat Jan 14 22:00:00 UTC 2006
9 //
10 //
11 
16 
28 
30  {
31  public:
32 
33  explicit TrackListMerger(const edm::ParameterSet& conf);
34 
35  virtual ~TrackListMerger();
36 
37  virtual void produce(edm::Event& e, const edm::EventSetup& c) override;
38 
39  private:
40  std::auto_ptr<reco::TrackCollection> outputTrks;
41  std::auto_ptr<reco::TrackExtraCollection> outputTrkExtras;
42  std::auto_ptr< TrackingRecHitCollection> outputTrkHits;
43  std::auto_ptr< std::vector<Trajectory> > outputTrajs;
44  std::auto_ptr< TrajTrackAssociationCollection > outputTTAss;
45  std::auto_ptr< TrajectorySeedCollection > outputSeeds;
46 
52 
55 
56  struct TkEDGetTokenss {
65  edm::EDGetTokenT<std::vector<Trajectory> > && traj_, edm::EDGetTokenT<TrajTrackAssociationCollection > && tass_,
67  tag(tag_), tk(tk_), traj(traj_), tass(tass_), tsel(tsel_), tmva(tmva_) {}
68  };
69  TkEDGetTokenss edTokens(const edm::InputTag &tag, const edm::InputTag &seltag, const edm::InputTag &mvatag) {
70  return TkEDGetTokenss(tag, consumes<reco::TrackCollection>(tag),
71  consumes<std::vector<Trajectory> >(tag), consumes<TrajTrackAssociationCollection >(tag),
72  consumes<edm::ValueMap<int> >(seltag), consumes<edm::ValueMap<float> >(mvatag));
73  }
75  return TkEDGetTokenss(tag, consumes<reco::TrackCollection>(tag),
76  consumes<std::vector<Trajectory> >(tag), consumes<TrajTrackAssociationCollection >(tag),
78  }
79  std::vector<TkEDGetTokenss> trackProducers_;
80 
82  double minPT_;
83  unsigned int minFound_;
84  float epsilon_;
85  float shareFrac_;
88  std::vector<double> indivShareFrac_;
89 
90  std::vector< std::vector< int> > listsToMerge_;
91  std::vector<bool> promoteQuality_;
92  std::vector<int> hasSelector_;
93  bool copyMVA_;
94 
99 
100  };
101 
102 
103 
104 #include <memory>
105 #include <string>
106 #include <iostream>
107 #include <cmath>
108 #include <vector>
109 
117 
119 
124 
125 //#include "DataFormats/TrackReco/src/classes.h"
126 
128 
129 #ifdef STAT_TSB
130 #include <x86intrin.h>
131 #endif
132 
133 namespace {
134 #ifdef STAT_TSB
135 inline volatile unsigned long long rdtsc() {
136  return __rdtsc();
137 }
138 
139  struct StatCount {
140  float maxDP=0.;
141  float maxDE=0.;
142  unsigned long long st;
143  long long totBegin=0;
144  long long totPre=0;
145  long long totEnd=0;
146  unsigned long long timeNo; // no-overlap
147  unsigned long long timeOv; // overlap
148  void begin(int tt) {
149  totBegin+=tt;
150  }
151  void start() { st=rdtsc(); }
152  void noOverlap() { timeNo += (rdtsc()-st);}
153  void overlap() { timeOv += (rdtsc()-st);}
154  void pre(int tt) { totPre+=tt;}
155  void end(int tt) { totEnd+=tt;}
156  void de(float d) { if (d>maxDE) maxDE=d;}
157  void dp(float d) { if (d>maxDP) maxDP=d;}
158 
159 
160  void print() const {
161  std::cout << "TrackListMerger stat\nBegin/Pre/End/maxDPhi/maxDEta/Overlap/NoOverlap "
162  << totBegin <<'/'<< totPre <<'/'<< totEnd <<'/'<< maxDP <<'/'<< maxDE
163  <<'/'<< timeOv/1000 <<'/'<< timeNo/1000
164  << std::endl;
165  }
166  StatCount() {}
167  ~StatCount() { print();}
168  };
169  StatCount statCount;
170 
171 #else
172  struct StatCount {
173  void begin(int){}
174  void pre(int){}
175  void end(int){}
176  void start(){}
177  void noOverlap(){}
178  void overlap(){}
179  void de(float){}
180  void dp(float){}
181 
182 
183  };
184  [[cms::thread_safe]] StatCount statCount;
185 #endif
186 
187 
188 }
189 
190 
191 
192 namespace {
193  edm::ProductID clusterProductB( const TrackingRecHit *hit){
194  return reinterpret_cast<const BaseTrackerRecHit *>(hit)->firstClusterRef().id();
195  }
196 }
197 
198 
200  copyExtras_ = conf.getUntrackedParameter<bool>("copyExtras", true);
201 
202  std::vector<edm::InputTag> trackProducerTags(conf.getParameter<std::vector<edm::InputTag> >("TrackProducers"));
203  //which of these do I need to turn into vectors?
204  maxNormalizedChisq_ = conf.getParameter<double>("MaxNormalizedChisq");
205  minPT_ = conf.getParameter<double>("MinPT");
206  minFound_ = (unsigned int)conf.getParameter<int>("MinFound");
207  epsilon_ = conf.getParameter<double>("Epsilon");
208  shareFrac_ = conf.getParameter<double>("ShareFrac");
209  allowFirstHitShare_ = conf.getParameter<bool>("allowFirstHitShare");
210  foundHitBonus_ = conf.getParameter<double>("FoundHitBonus");
211  lostHitPenalty_ = conf.getParameter<double>("LostHitPenalty");
212  indivShareFrac_=conf.getParameter<std::vector<double> >("indivShareFrac");
213  std::string qualityStr = conf.getParameter<std::string>("newQuality");
214 
215  if (qualityStr != "") {
217  }
218  else
220 
221  use_sharesInput_ = true;
222  if ( epsilon_ > 0.0 )use_sharesInput_ = false;
223 
224  edm::VParameterSet setsToMerge=conf.getParameter<edm::VParameterSet>("setsToMerge");
225 
226  for ( unsigned int i=0; i<setsToMerge.size(); i++) {
227  listsToMerge_.push_back(setsToMerge[i].getParameter<std::vector< int> >("tLists"));
228  promoteQuality_.push_back(setsToMerge[i].getParameter<bool>("pQual"));
229  }
230  hasSelector_ = conf.getParameter<std::vector<int> >("hasSelector");
231  copyMVA_ = conf.getParameter<bool>("copyMVA");
232 
233  std::vector<edm::InputTag> selectors(conf.getParameter<std::vector<edm::InputTag> >("selectedTrackQuals"));
234  std::vector<edm::InputTag> mvaStores;
235  if(conf.exists("mvaValueTags")){
236  mvaStores = conf.getParameter<std::vector<edm::InputTag> >("mvaValueTags");
237  }else{
238  for (int i = 0; i < (int)selectors.size(); i++){
239  edm::InputTag ntag(selectors[i].label(),"MVAVals");
240  mvaStores.push_back(ntag);
241  }
242  }
243  unsigned int numTrkColl = trackProducerTags.size();
244  if (numTrkColl != hasSelector_.size() || numTrkColl != selectors.size()) {
245  throw cms::Exception("Inconsistent size") << "need same number of track collections and selectors";
246  }
247  if (numTrkColl != hasSelector_.size() || numTrkColl != mvaStores.size()) {
248  throw cms::Exception("Inconsistent size") << "need same number of track collections and MVA stores";
249  }
250  for (unsigned int i = indivShareFrac_.size(); i < numTrkColl; i++) {
251  // edm::LogWarning("TrackListMerger") << "No indivShareFrac for " << trackProducersTags <<". Using default value of 1";
252  indivShareFrac_.push_back(1.0);
253  }
254 
255  trkQualMod_=conf.getParameter<bool>("writeOnlyTrkQuals");
256  if ( trkQualMod_) {
257  bool ok=true;
258  for ( unsigned int i=1; i<numTrkColl; i++) {
259  if (!(trackProducerTags[i]==trackProducerTags[0])) ok=false;
260  }
261  if ( !ok) {
262  throw cms::Exception("Bad input") << "to use writeOnlyTrkQuals=True all input InputTags must be the same";
263  }
264  produces<edm::ValueMap<int> >();
265  }
266  else{
267  produces<reco::TrackCollection>();
268 
269  makeReKeyedSeeds_ = conf.getUntrackedParameter<bool>("makeReKeyedSeeds",false);
270  if (makeReKeyedSeeds_){
271  copyExtras_=true;
272  produces<TrajectorySeedCollection>();
273  }
274 
275  if (copyExtras_) {
276  produces<reco::TrackExtraCollection>();
277  produces<TrackingRecHitCollection>();
278  }
279  produces< std::vector<Trajectory> >();
280  produces< TrajTrackAssociationCollection >();
281  }
282  produces<edm::ValueMap<float> >("MVAVals");
283 
284  // Do all the consumes
285  trackProducers_.resize(numTrkColl);
286  for (unsigned int i = 0; i < numTrkColl; ++i) {
287  trackProducers_[i] = hasSelector_[i]>0 ? edTokens(trackProducerTags[i], selectors[i], mvaStores[i]) : edTokens(trackProducerTags[i], mvaStores[i]);
288  }
289 }
290 
291 
292 // Virtual destructor needed.
294 
295 // Functions that gets called by framework every event
297 {
298  // extract tracker geometry
299  //
300  //edm::ESHandle<TrackerGeometry> theG;
301  //es.get<TrackerDigiGeometryRecord>().get(theG);
302 
303  // using namespace reco;
304 
305  // get Inputs
306  // if 1 input list doesn't exist, make an empty list, issue a warning, and continue
307  // this allows TrackListMerger to be used as a cleaner only if handed just one list
308  // if both input lists don't exist, will issue 2 warnings and generate an empty output collection
309  //
310  static const reco::TrackCollection s_empty;
311 
312  std::vector<const reco::TrackCollection *> trackColls;
313  std::vector<edm::Handle<reco::TrackCollection> > trackHandles(trackProducers_.size());
314  for ( unsigned int i=0; i<trackProducers_.size(); i++) {
315  trackColls.push_back(0);
316  //edm::Handle<reco::TrackCollection> trackColl;
317  e.getByToken(trackProducers_[i].tk, trackHandles[i]);
318  if (trackHandles[i].isValid()) {
319  trackColls[i]= trackHandles[i].product();
320  } else {
321  edm::LogWarning("TrackListMerger") << "TrackCollection " << trackProducers_[i].tag <<" not found";
322  trackColls[i]=&s_empty;
323  }
324  }
325 
326  unsigned int collsSize =trackColls.size();
327  unsigned int rSize=0;
328  unsigned int trackCollSizes[collsSize];
329  unsigned int trackCollFirsts[collsSize];
330  for (unsigned int i=0; i!=collsSize; i++) {
331  trackCollSizes[i]=trackColls[i]->size();
332  trackCollFirsts[i]=rSize;
333  rSize+=trackCollSizes[i];
334  }
335 
336  statCount.begin(rSize);
337 
338  //
339  // quality cuts first
340  //
341  int i=-1;
342 
343  int selected[rSize];
344  int indexG[rSize];
345  bool trkUpdated[rSize];
346  int trackCollNum[rSize];
347  int trackQuals[rSize];
348  float trackMVAs[rSize];
349  for (unsigned int j=0; j<rSize;j++) {
350  indexG[j]=-1; selected[j]=1; trkUpdated[j]=false; trackCollNum[j]=0; trackQuals[j]=0;trackMVAs[j] = -99.0;
351  }
352 
353  int ngood=0;
354  for (unsigned int j=0; j!= collsSize; j++) {
355  const reco::TrackCollection *tC1=trackColls[j];
356 
357  edm::Handle<edm::ValueMap<int> > trackSelColl;
358  edm::Handle<edm::ValueMap<float> > trackMVAStore;
359  if ( copyMVA_ )
360  e.getByToken(trackProducers_[j].tmva, trackMVAStore);
361  if ( hasSelector_[j]>0 ){
362  e.getByToken(trackProducers_[j].tsel, trackSelColl);
363  }
364 
365  if ( 0<tC1->size() ){
366  unsigned int iC=0;
367  for (reco::TrackCollection::const_iterator track=tC1->begin(); track!=tC1->end(); track++){
368  i++;
369  trackCollNum[i]=j;
370  trackQuals[i]=track->qualityMask();
371 
372  reco::TrackRef trkRef=reco::TrackRef(trackHandles[j],iC);
373  if ( copyMVA_ )
374  if( (*trackMVAStore).contains(trkRef.id()) ) trackMVAs[i] = (*trackMVAStore)[trkRef];
375  if ( hasSelector_[j]>0 ) {
376  int qual=(*trackSelColl)[trkRef];
377  if ( qual < 0 ) {
378  selected[i]=0;
379  iC++;
380  continue;
381  }else{
382  trackQuals[i]=qual;
383  }
384  }
385  iC++;
386  selected[i]=trackQuals[i]+10;//10 is magic number used throughout...
387  if ((short unsigned)track->ndof() < 1){
388  selected[i]=0;
389  continue;
390  }
391  if (track->normalizedChi2() > maxNormalizedChisq_){
392  selected[i]=0;
393  continue;
394  }
395  if (track->found() < minFound_){
396  selected[i]=0;
397  continue;
398  }
399  if (track->pt() < minPT_){
400  selected[i]=0;
401  continue;
402  }
403  // good!
404  indexG[i] = ngood++;
405  //if ( beVerb) std::cout << "inverb " << track->pt() << " " << selected[i] << std::endl;
406  }//end loop over tracks
407  }//end more than 0 track
408  } // loop over trackcolls
409 
410 
411  statCount.pre(ngood);
412 
413  //cache the id and rechits of valid hits
414  typedef std::pair<unsigned int, const TrackingRecHit*> IHit;
415  std::vector<IHit> rh1[ngood]; // an array of vectors!
416  //const TrackingRecHit* fh1[ngood]; // first hit...
417  uint8_t algo[ngood];
418  float score[ngood];
419 
420 
421  for ( unsigned int j=0; j<rSize; j++) {
422  if (selected[j]==0) continue;
423  int i = indexG[j];
424  assert(i>=0);
425  unsigned int collNum=trackCollNum[j];
426  unsigned int trackNum=j-trackCollFirsts[collNum];
427  const reco::Track *track=&((trackColls[collNum])->at(trackNum));
428 
429  algo[i]=track->algo();
430  int validHits=track->numberOfValidHits();
431  int validPixelHits=track->hitPattern().numberOfValidPixelHits();
432  int lostHits=track->numberOfLostHits();
433  score[i] = foundHitBonus_*validPixelHits+foundHitBonus_*validHits - lostHitPenalty_*lostHits - track->chi2();
434 
435 
436  rh1[i].reserve(validHits) ;
437  auto compById = [](IHit const & h1, IHit const & h2) {return h1.first < h2.first;};
438  // fh1[i] = &(**track->recHitsBegin());
439  for (trackingRecHit_iterator it = track->recHitsBegin(); it != track->recHitsEnd(); ++it) {
440  const TrackingRecHit* hit = &(**it);
441  unsigned int id = hit->rawId() ;
442  if(hit->geographicalId().subdetId()>2) id &= (~3); // mask mono/stereo in strips...
443  if likely(hit->isValid()) { rh1[i].emplace_back(id,hit); std::push_heap(rh1[i].begin(),rh1[i].end(),compById); }
444  }
445  std::sort_heap(rh1[i].begin(),rh1[i].end(),compById);
446  }
447 
448  //DL here
449  if likely(ngood>1 && collsSize>1)
450  for ( unsigned int ltm=0; ltm<listsToMerge_.size(); ltm++) {
451  int saveSelected[rSize];
452  bool notActive[collsSize];
453  for (unsigned int cn=0;cn!=collsSize;++cn)
454  notActive[cn]= find(listsToMerge_[ltm].begin(),listsToMerge_[ltm].end(),cn)==listsToMerge_[ltm].end();
455 
456  for ( unsigned int i=0; i<rSize; i++) saveSelected[i]=selected[i];
457 
458  //DL protect against 0 tracks?
459  for ( unsigned int i=0; i<rSize-1; i++) {
460  if (selected[i]==0) continue;
461  unsigned int collNum=trackCollNum[i];
462 
463  //check that this track is in one of the lists for this iteration
464  if (notActive[collNum]) continue;
465 
466  int k1 = indexG[i];
467  unsigned int nh1=rh1[k1].size();
468  int qualityMaskT1 = trackQuals[i];
469 
470  int nhit1 = nh1; // validHits[k1];
471  float score1 = score[k1];
472 
473  // start at next collection
474  for ( unsigned int j=i+1; j<rSize; j++) {
475  if (selected[j]==0) continue;
476  unsigned int collNum2=trackCollNum[j];
477  if ( (collNum == collNum2) && indivShareFrac_[collNum] > 0.99) continue;
478  //check that this track is in one of the lists for this iteration
479  if (notActive[collNum2]) continue;
480 
481  int k2 = indexG[j];
482 
483  int newQualityMask = -9; //avoid resetting quality mask if not desired 10+ -9 =1
484  if (promoteQuality_[ltm]) {
485  int maskT1= saveSelected[i]>1? saveSelected[i]-10 : qualityMaskT1;
486  int maskT2= saveSelected[j]>1? saveSelected[j]-10 : trackQuals[j];
487  newQualityMask =(maskT1 | maskT2); // take OR of trackQuality
488  }
489  unsigned int nh2=rh1[k2].size();
490  int nhit2 = nh2;
491 
492 
493  auto share = use_sharesInput_ ?
494  [](const TrackingRecHit* it,const TrackingRecHit* jt, float)->bool { return it->sharesInput(jt,TrackingRecHit::some); } :
495  [](const TrackingRecHit* it,const TrackingRecHit* jt, float eps)->bool {
496  float delta = std::abs ( it->localPosition().x()-jt->localPosition().x() );
497  return (it->geographicalId()==jt->geographicalId())&&(delta<eps);
498  };
499 
500  statCount.start();
501 
502  //loop over rechits
503  int noverlap=0;
504  int firstoverlap=0;
505  // check first hit (should use REAL first hit?)
506  if unlikely(allowFirstHitShare_ && rh1[k1][0].first==rh1[k2][0].first ) {
507  const TrackingRecHit* it = rh1[k1][0].second;
508  const TrackingRecHit* jt = rh1[k2][0].second;
509  if (share(it,jt,epsilon_)) firstoverlap=1;
510  }
511 
512 
513  // exploit sorting
514  unsigned int jh=0;
515  unsigned int ih=0;
516  while (ih!=nh1 && jh!=nh2) {
517  // break if not enough to go...
518  // if ( nprecut-noverlap+firstoverlap > int(nh1-ih)) break;
519  // if ( nprecut-noverlap+firstoverlap > int(nh2-jh)) break;
520  auto const id1 = rh1[k1][ih].first;
521  auto const id2 = rh1[k2][jh].first;
522  if (id1<id2) ++ih;
523  else if (id2<id1) ++jh;
524  else {
525  // in case of split-hit do full conbinatorics
526  auto li=ih; while( (++li)!=nh1 && id1 == rh1[k1][li].first);
527  auto lj=jh; while( (++lj)!=nh2 && id2 == rh1[k2][lj].first);
528  for (auto ii=ih; ii!=li; ++ii)
529  for (auto jj=jh; jj!=lj; ++jj) {
530  const TrackingRecHit* it = rh1[k1][ii].second;
531  const TrackingRecHit* jt = rh1[k2][jj].second;
532  if (share(it,jt,epsilon_)) noverlap++;
533  }
534  jh=lj; ih=li;
535  } // equal ids
536 
537  } //loop over ih & jh
538 
539  bool dupfound = (collNum != collNum2) ? (noverlap-firstoverlap) > (std::min(nhit1,nhit2)-firstoverlap)*shareFrac_ :
540  (noverlap-firstoverlap) > (std::min(nhit1,nhit2)-firstoverlap)*indivShareFrac_[collNum];
541 
542  if ( dupfound ) {
543  float score2 = score[k2];
544  constexpr float almostSame = 0.01f; // difference rather than ratio due to possible negative values for score
545  if ( score1 - score2 > almostSame ) {
546  selected[j]=0;
547  selected[i]=10+newQualityMask; // add 10 to avoid the case where mask = 1
548  trkUpdated[i]=true;
549  } else if ( score2 - score1 > almostSame ) {
550  selected[i]=0;
551  selected[j]=10+newQualityMask; // add 10 to avoid the case where mask = 1
552  trkUpdated[j]=true;
553  }else{
554  // If tracks from both iterations are virtually identical, choose the one with the best quality or with lower algo
557  //same quality, pick earlier algo
558  if (algo[k1] <= algo[k2]) {
559  selected[j]=0;
560  selected[i]=10+newQualityMask; // add 10 to avoid the case where mask = 1
561  trkUpdated[i]=true;
562  } else {
563  selected[i]=0;
564  selected[j]=10+newQualityMask; // add 10 to avoid the case where mask = 1
565  trkUpdated[j]=true;
566  }
567  } else if ((trackQuals[j] & (1<<reco::TrackBase::loose|1<<reco::TrackBase::tight|1<<reco::TrackBase::highPurity) ) <
569  selected[j]=0;
570  selected[i]=10+newQualityMask; // add 10 to avoid the case where mask = 1
571  trkUpdated[i]=true;
572  }else{
573  selected[i]=0;
574  selected[j]=10+newQualityMask; // add 10 to avoid the case where mask = 1
575  trkUpdated[j]=true;
576  }
577  }//end fi < fj
578  statCount.overlap();
579  /*
580  if (at0[k1]&&at0[k2]) {
581  statCount.dp(dphi);
582  if (dz<1.f) statCount.de(deta);
583  }
584  */
585  }//end got a duplicate
586  else {
587  statCount.noOverlap();
588  }
589  //stop if the ith track is now unselected
590  if (selected[i]==0) break;
591  }//end track2 loop
592  }//end track loop
593  } //end loop over track list sets
594 
595 
596 
597  std::auto_ptr<edm::ValueMap<float> > vmMVA = std::auto_ptr<edm::ValueMap<float> >(new edm::ValueMap<float>);
598  edm::ValueMap<float>::Filler fillerMVA(*vmMVA);
599 
600 
601 
602  // special case - if just doing the trkquals
603  if (trkQualMod_) {
604  unsigned int tSize=trackColls[0]->size();
605  std::auto_ptr<edm::ValueMap<int> > vm = std::auto_ptr<edm::ValueMap<int> >(new edm::ValueMap<int>);
606  edm::ValueMap<int>::Filler filler(*vm);
607 
608  std::vector<int> finalQuals(tSize,-1); //default is unselected
609  for ( unsigned int i=0; i<rSize; i++) {
610  unsigned int tNum=i%tSize;
611 
612  if (selected[i]>1 ) {
613  finalQuals[tNum]=selected[i]-10;
614  if (trkUpdated[i])
615  finalQuals[tNum]=(finalQuals[tNum] | (1<<qualityToSet_));
616  }
617  if ( selected[i]==1 )
618  finalQuals[tNum]=trackQuals[i];
619  }
620 
621  filler.insert(trackHandles[0], finalQuals.begin(),finalQuals.end());
622  filler.fill();
623 
624  e.put(vm);
625 
626  std::vector<float> mvaVec(tSize,-99);
627 
628  for(unsigned int i = 0; i < rSize; i++){
629  unsigned int tNum = i % tSize;
630  mvaVec[tNum] = trackMVAs[tNum];
631  }
632 
633  fillerMVA.insert(trackHandles[0],mvaVec.begin(),mvaVec.end());
634  fillerMVA.fill();
635  if ( copyMVA_)
636  e.put(vmMVA,"MVAVals");
637  return;
638  }
639 
640 
641  //
642  // output selected tracks - if any
643  //
644 
645  reco::TrackRef trackRefs[rSize];
646  edm::RefToBase<TrajectorySeed> seedsRefs[rSize];
647 
648  unsigned int nToWrite=0;
649  for ( unsigned int i=0; i<rSize; i++)
650  if (selected[i]!=0) nToWrite++;
651 
652  std::vector<float> mvaVec;
653 
654  outputTrks = std::auto_ptr<reco::TrackCollection>(new reco::TrackCollection);
655  outputTrks->reserve(nToWrite);
657 
658  if (copyExtras_) {
659  outputTrkExtras = std::auto_ptr<reco::TrackExtraCollection>(new reco::TrackExtraCollection);
660  outputTrkExtras->reserve(nToWrite);
662  outputTrkHits = std::auto_ptr<TrackingRecHitCollection>(new TrackingRecHitCollection);
663  outputTrkHits->reserve(nToWrite*25);
665  if (makeReKeyedSeeds_){
666  outputSeeds = std::auto_ptr<TrajectorySeedCollection>(new TrajectorySeedCollection);
667  outputSeeds->reserve(nToWrite);
669  }
670  }
671 
672 
673  outputTrajs = std::auto_ptr< std::vector<Trajectory> >(new std::vector<Trajectory>());
674  outputTrajs->reserve(rSize);
675  outputTTAss = std::auto_ptr< TrajTrackAssociationCollection >(new TrajTrackAssociationCollection());
676 
677 
678 
679  for ( unsigned int i=0; i<rSize; i++) {
680  if (selected[i]==0) {
681  trackRefs[i]=reco::TrackRef();
682  continue;
683  }
684 
685  unsigned int collNum=trackCollNum[i];
686  unsigned int trackNum=i-trackCollFirsts[collNum];
687  const reco::Track *track=&((trackColls[collNum])->at(trackNum));
688  outputTrks->push_back( reco::Track( *track ) );
689  mvaVec.push_back(trackMVAs[i]);
690  if (selected[i]>1 ) {
691  outputTrks->back().setQualityMask(selected[i]-10);
692  if (trkUpdated[i])
693  outputTrks->back().setQuality(qualityToSet_);
694  }
695  //might duplicate things, but doesnt hurt
696  if ( selected[i]==1 )
697  outputTrks->back().setQualityMask(trackQuals[i]);
698 
699  // if ( beVerb ) std::cout << "selected " << outputTrks->back().pt() << " " << outputTrks->back().qualityMask() << " " << selected[i] << std::endl;
700 
701  //fill the TrackCollection
702  if (copyExtras_) {
703  edm::RefToBase<TrajectorySeed> origSeedRef = track->seedRef();
704  //creating a seed with rekeyed clusters if required
705  if (makeReKeyedSeeds_){
706  bool doRekeyOnThisSeed=false;
707 
708  edm::InputTag clusterRemovalInfos("");
709  //grab on of the hits of the seed
710  if (origSeedRef->nHits()!=0){
711  TrajectorySeed::const_iterator firstHit=origSeedRef->recHits().first;
712  const TrackingRecHit *hit = &*firstHit;
713  if (firstHit->isValid()){
714  edm::ProductID pID=clusterProductB(hit);
715  // the cluster collection either produced a removalInfo or mot
716  //get the clusterremoval info from the provenance: will rekey if this is found
718  edm::Provenance prov=e.getProvenance(pID);
719  clusterRemovalInfos=edm::InputTag(prov.moduleLabel(),
720  prov.productInstanceName(),
721  prov.processName());
722  doRekeyOnThisSeed=e.getByLabel(clusterRemovalInfos,CRIh);
723  }//valid hit
724  }//nhit!=0
725 
726  if (doRekeyOnThisSeed && !(clusterRemovalInfos==edm::InputTag(""))) {
727  ClusterRemovalRefSetter refSetter(e,clusterRemovalInfos);
728  TrajectorySeed::recHitContainer newRecHitContainer;
729  newRecHitContainer.reserve(origSeedRef->nHits());
730  TrajectorySeed::const_iterator iH=origSeedRef->recHits().first;
731  TrajectorySeed::const_iterator iH_end=origSeedRef->recHits().second;
732  for (;iH!=iH_end;++iH){
733  newRecHitContainer.push_back(*iH);
734  refSetter.reKey(&newRecHitContainer.back());
735  }
736  outputSeeds->push_back( TrajectorySeed( origSeedRef->startingState(),
737  newRecHitContainer,
738  origSeedRef->direction()));
739  }
740  //doRekeyOnThisSeed=true
741  else{
742  //just copy the one we had before
743  outputSeeds->push_back( TrajectorySeed(*origSeedRef));
744  }
746  origSeedRef=edm::RefToBase<TrajectorySeed>( pureRef);
747  }//creating a new seed and rekeying it rechit clusters.
748 
749  // Fill TrackExtra collection
750  outputTrkExtras->push_back( reco::TrackExtra(
751  track->outerPosition(), track->outerMomentum(), track->outerOk(),
752  track->innerPosition(), track->innerMomentum(), track->innerOk(),
753  track->outerStateCovariance(), track->outerDetId(),
754  track->innerStateCovariance(), track->innerDetId(),
755  track->seedDirection(), origSeedRef ) );
756  seedsRefs[i]=origSeedRef;
757  outputTrks->back().setExtra( reco::TrackExtraRef( refTrkExtras, outputTrkExtras->size() - 1) );
758  reco::TrackExtra & tx = outputTrkExtras->back();
759  tx.setResiduals(track->residuals());
760 
761  // fill TrackingRecHits
762  unsigned nh1=track->recHitsSize();
763  for ( unsigned ih=0; ih<nh1; ++ih ) {
764  //const TrackingRecHit*hit=&((*(track->recHit(ih))));
765  outputTrkHits->push_back( track->recHit(ih)->clone() );
766  tx.add( TrackingRecHitRef( refTrkHits, outputTrkHits->size() - 1) );
767  }
768  }
769  trackRefs[i] = reco::TrackRef(refTrks, outputTrks->size() - 1);
770 
771 
772  }//end faux loop over tracks
773 
774  //Fill the trajectories, etc. for 1st collection
775  refTrajs = e.getRefBeforePut< std::vector<Trajectory> >();
776 
777  for (unsigned int ti=0; ti<trackColls.size(); ti++) {
780  e.getByToken(trackProducers_[ti].traj, hTraj1);
781  e.getByToken(trackProducers_[ti].tass, hTTAss1);
782 
783  if (hTraj1.failedToGet() || hTTAss1.failedToGet()) continue;
784 
785  for (size_t i = 0, n = hTraj1->size(); i < n; ++i) {
786  edm::Ref< std::vector<Trajectory> > trajRef(hTraj1, i);
787  TrajTrackAssociationCollection::const_iterator match = hTTAss1->find(trajRef);
788  if (match != hTTAss1->end()) {
789  const edm::Ref<reco::TrackCollection> &trkRef = match->val;
790  short oldKey = trackCollFirsts[ti]+static_cast<short>(trkRef.key());
791  if (trackRefs[oldKey].isNonnull()) {
792  outputTrajs->push_back( *trajRef );
793  //if making extras and the seeds at the same time, change the seed ref on the trajectory
795  outputTrajs->back().setSeedRef( seedsRefs[oldKey] );
796  outputTTAss->insert ( edm::Ref< std::vector<Trajectory> >(refTrajs, outputTrajs->size() - 1),
797  trackRefs[oldKey] );
798  }
799  }
800  }
801  }
802 
803  statCount.end(outputTrks->size());
804 
805  edm::ProductID nPID = refTrks.id();
807  fillerMVA.insert(outHandle,mvaVec.begin(),mvaVec.end());
808  fillerMVA.fill();
809 
810  e.put(outputTrks);
811  if ( copyMVA_ )
812  e.put(vmMVA,"MVAVals");
813  if (copyExtras_) {
814  e.put(outputTrkExtras);
815  e.put(outputTrkHits);
816  if (makeReKeyedSeeds_)
817  e.put(outputSeeds);
818  }
819  e.put(outputTrajs);
820  e.put(outputTTAss);
821  return;
822 
823 }//end produce
824 
825 
828 
830 
PropagationDirection direction() const
std::auto_ptr< reco::TrackCollection > outputTrks
dbl * delta
Definition: mlp_gen.cc:36
T getParameter(std::string const &) const
#define dso_hidden
T getUntrackedParameter(std::string const &, T const &) const
int i
Definition: DBlmapReader.cc:9
TrackListMerger(const edm::ParameterSet &conf)
reference back()
Definition: OwnVector.h:327
tuple start
Check for commandline option errors.
Definition: dqm_diff.py:58
std::vector< bool > promoteQuality_
reco::TrackRefProd refTrks
bool getByToken(EDGetToken token, Handle< PROD > &result) const
Definition: Event.h:446
edm::EDGetTokenT< TrajTrackAssociationCollection > tass
edm::RefProd< std::vector< Trajectory > > refTrajs
TkEDGetTokenss edTokens(const edm::InputTag &tag, const edm::InputTag &mvatag)
virtual void produce(Event &, EventSetup const &)=0
TrackQuality
track quality
Definition: TrackBase.h:113
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:17
size_t recHitsSize() const
Get number of RecHits. (Warning, this includes invalid hits, which are not physical hits)...
Definition: Track.h:119
std::vector< ParameterSet > VParameterSet
Definition: ParameterSet.h:33
std::string print(const Track &, edm::Verbosity=edm::Concise)
Track print utility.
Definition: print.cc:10
void reKey(TrackingRecHit *hit) const
reco::TrackBase::TrackQuality qualityToSet_
std::vector< Track > TrackCollection
collection of Tracks
Definition: TrackFwd.h:13
bool exists(std::string const &parameterName) const
checks if a parameter exists
bool innerOk() const
return true if the innermost hit is valid
Definition: Track.h:50
unsigned short numberOfLostHits() const
number of cases where track crossed a layer without getting a hit.
Definition: TrackBase.h:802
key_type key() const
Accessor for product key.
Definition: Ref.h:266
int ii
Definition: cuy.py:588
std::string const & processName() const
Definition: Provenance.h:62
void find(edm::Handle< EcalRecHitCollection > &hits, DetId thisDet, std::vector< EcalRecHitCollection::const_iterator > &hit, bool debug=false)
Definition: FindCaloHit.cc:7
const math::XYZPoint & outerPosition() const
position of the outermost hit
Definition: Track.h:65
edm::RefProd< TrajectorySeedCollection > refTrajSeeds
#define constexpr
ProductID id() const
Accessor for product ID.
Definition: Ref.h:256
TrackingRecHitRefProd refTrkHits
virtual void produce(edm::Event &e, const edm::EventSetup &c) override
std::vector< double > indivShareFrac_
bool overlap(const reco::Muon &muon1, const reco::Muon &muon2, double pullX=1.0, double pullY=1.0, bool checkAdjacentChambers=false)
const math::XYZPoint & innerPosition() const
position of the innermost hit
Definition: Track.h:55
TrackAlgorithm algo() const
Definition: TrackBase.h:403
#define unlikely(x)
#define likely(x)
void push_back(D *&d)
Definition: OwnVector.h:280
std::auto_ptr< TrackingRecHitCollection > outputTrkHits
edm::EDGetTokenT< edm::ValueMap< int > > tsel
std::vector< TrajectorySeed > TrajectorySeedCollection
double chi2() const
chi-squared of the fit
Definition: TrackBase.h:525
OrphanHandle< PROD > put(std::auto_ptr< PROD > product)
Put a new product.
Definition: Event.h:113
recHitContainer::const_iterator const_iterator
edm::Ref< TrackingRecHitCollection > TrackingRecHitRef
persistent reference to a TrackingRecHit
CovarianceMatrix outerStateCovariance() const
outermost trajectory state curvilinear errors
Definition: Track.h:75
edm::EDGetTokenT< reco::TrackCollection > tk
unsigned int outerDetId() const
DetId of the detector on which surface the outermost state is located.
Definition: Track.h:94
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
int j
Definition: DBlmapReader.cc:9
TkEDGetTokenss edTokens(const edm::InputTag &tag, const edm::InputTag &seltag, const edm::InputTag &mvatag)
unsigned short numberOfValidHits() const
number of valid hits found
Definition: TrackBase.h:796
#define end
Definition: vmac.h:37
std::auto_ptr< TrajTrackAssociationCollection > outputTTAss
trackingRecHit_iterator recHitsBegin() const
Iterator to first hit on the track.
Definition: Track.h:104
virtual ~TrackListMerger()
bool first
Definition: L1TdeRCT.cc:75
edm::AssociationMap< edm::OneToOne< std::vector< Trajectory >, reco::TrackCollection, unsigned short > > TrajTrackAssociationCollection
int subdetId() const
get the contents of the subdetector field (not cast into any detector&#39;s numbering enum) ...
Definition: DetId.h:37
std::auto_ptr< std::vector< Trajectory > > outputTrajs
std::auto_ptr< TrajectorySeedCollection > outputSeeds
bool getByLabel(InputTag const &tag, Handle< PROD > &result) const
Definition: Event.h:402
RefProd< PROD > getRefBeforePut()
Definition: Event.h:133
tuple conf
Definition: dbtoconf.py:185
std::vector< int > hasSelector_
bool failedToGet() const
Definition: HandleBase.h:80
std::vector< TrackExtra > TrackExtraCollection
collection of TrackExtra objects
Definition: TrackExtraFwd.h:11
double maxNormalizedChisq_
edm::OwnVector< TrackingRecHit > TrackingRecHitCollection
collection of TrackingRecHits
edm::RefToBase< TrajectorySeed > seedRef() const
Definition: Track.h:213
PTrajectoryStateOnDet const & startingState() const
static TrackQuality qualityByName(const std::string &name)
Definition: TrackBase.cc:102
const math::XYZVector & outerMomentum() const
momentum vector at the outermost hit position
Definition: Track.h:70
bool outerOk() const
return true if the outermost hit is valid
Definition: Track.h:45
const HitPattern & hitPattern() const
Access the hit pattern, indicating in which Tracker layers the track has hits.
Definition: TrackBase.h:364
edm::EDGetTokenT< edm::ValueMap< float > > tmva
edm::Ref< TrackCollection > TrackRef
persistent reference to a Track
Definition: TrackFwd.h:19
CovarianceMatrix innerStateCovariance() const
innermost trajectory state curvilinear errors
Definition: Track.h:80
range recHits() const
bool isValid() const
void add(const TrackingRecHitRef &r)
add a reference to a RecHit
std::string const & moduleLabel() const
Definition: Provenance.h:61
std::auto_ptr< reco::TrackExtraCollection > outputTrkExtras
unsigned int nHits() const
#define begin
Definition: vmac.h:30
unsigned int minFound_
std::vector< TkEDGetTokenss > trackProducers_
const math::XYZVector & innerMomentum() const
momentum vector at the innermost hit position
Definition: Track.h:60
TrackingRecHitRef recHit(size_t i) const
Get i-th hit on the track.
Definition: Track.h:114
ProductID id() const
Accessor for product ID.
Definition: RefProd.h:140
const TrackResiduals & residuals() const
Definition: Track.h:220
int numberOfValidPixelHits() const
Definition: HitPattern.h:752
edm::EDGetTokenT< std::vector< Trajectory > > traj
PropagationDirection seedDirection() const
direction of how the hits were sorted in the original seed
Definition: Track.h:204
tuple cout
Definition: gather_cfg.py:121
std::pair< typename Association::data_type::first_type, double > match(Reference key, Association association, bool bestMatchByMaxValue)
Generic matching function.
Definition: Utils.h:6
TkEDGetTokenss(const edm::InputTag &tag_, edm::EDGetTokenT< reco::TrackCollection > &&tk_, edm::EDGetTokenT< std::vector< Trajectory > > &&traj_, edm::EDGetTokenT< TrajTrackAssociationCollection > &&tass_, edm::EDGetTokenT< edm::ValueMap< int > > &&tsel_, edm::EDGetTokenT< edm::ValueMap< float > > &&tmva_)
DetId geographicalId() const
std::string const & productInstanceName() const
Definition: Provenance.h:63
std::vector< std::vector< int > > listsToMerge_
Provenance getProvenance(BranchID const &theID) const
Definition: Event.cc:73
T x() const
Definition: PV3DBase.h:62
virtual LocalPoint localPosition() const =0
unsigned int innerDetId() const
DetId of the detector on which surface the innermost state is located.
Definition: Track.h:99
id_type rawId() const
void reserve(size_t)
Definition: OwnVector.h:274
reco::TrackExtraRefProd refTrkExtras
list at
Definition: asciidump.py:428
trackingRecHit_iterator recHitsEnd() const
Iterator to last hit on the track.
Definition: Track.h:109
void setResiduals(const TrackResiduals &r)
set the residuals
Definition: TrackExtra.h:178