Inheritance diagram for TrackListMerger:

Classes
struct	TkEDGetTokenss

Public Member Functions
virtual void	produce (edm::Event &e, const edm::EventSetup &c) override

	TrackListMerger (const edm::ParameterSet &conf)

virtual	~TrackListMerger ()

Public Member Functions inherited from edm::stream::EDProducer<>
	EDProducer ()=default

Public Member Functions inherited from edm::stream::EDProducerBase
	EDProducerBase ()

ModuleDescription const &	moduleDescription () const

virtual	~EDProducerBase ()

Public Member Functions inherited from edm::ProducerBase
	ProducerBase ()

void	registerProducts (ProducerBase , ProductRegistry , ModuleDescription const &)

std::function< void(BranchDescription const &)>	registrationCallback () const
	used by the fwk to register list of products More...

virtual	~ProducerBase ()

Public Member Functions inherited from edm::EDConsumerBase
	EDConsumerBase ()

ProductHolderIndexAndSkipBit	indexFrom (EDGetToken, BranchType, TypeID const &) const

void	itemsMayGet (BranchType, std::vector< ProductHolderIndexAndSkipBit > &) const

void	itemsToGet (BranchType, std::vector< ProductHolderIndexAndSkipBit > &) const

std::vector < ProductHolderIndexAndSkipBit > const &	itemsToGetFromEvent () const

void	labelsForToken (EDGetToken iToken, Labels &oLabels) const

void	modulesDependentUpon (const std::string &iProcessName, std::vector< const char * > &oModuleLabels) const

bool	registeredToConsume (ProductHolderIndex, bool, BranchType) const

bool	registeredToConsumeMany (TypeID const &, BranchType) const

void	updateLookup (BranchType iBranchType, ProductHolderIndexHelper const &)

virtual	~EDConsumerBase ()

Private Member Functions
TkEDGetTokenss	edTokens (const edm::InputTag &tag, const edm::InputTag &seltag, const edm::InputTag &mvatag)

TkEDGetTokenss	edTokens (const edm::InputTag &tag, const edm::InputTag &mvatag)

Private Attributes
bool	allowFirstHitShare_

bool	copyExtras_

bool	copyMVA_

float	epsilon_

float	foundHitBonus_

std::vector< int >	hasSelector_

std::vector< double >	indivShareFrac_

std::vector< std::vector< int > >	listsToMerge_

float	lostHitPenalty_

bool	makeReKeyedSeeds_

double	maxNormalizedChisq_

unsigned int	minFound_

double	minPT_

std::auto_ptr < TrajectorySeedCollection >	outputSeeds

std::auto_ptr< std::vector < Trajectory > >	outputTrajs

std::auto_ptr < reco::TrackExtraCollection >	outputTrkExtras

std::auto_ptr < TrackingRecHitCollection >	outputTrkHits

std::auto_ptr < reco::TrackCollection >	outputTrks

std::auto_ptr < TrajTrackAssociationCollection >	outputTTAss

std::vector< bool >	promoteQuality_

reco::TrackBase::TrackQuality	qualityToSet_

edm::RefProd< std::vector < Trajectory > >	refTrajs

edm::RefProd < TrajectorySeedCollection >	refTrajSeeds

reco::TrackExtraRefProd	refTrkExtras

TrackingRecHitRefProd	refTrkHits

reco::TrackRefProd	refTrks

float	shareFrac_

std::vector< TkEDGetTokenss >	trackProducers_

bool	trkQualMod_

bool	use_sharesInput_

Additional Inherited Members
Public Types inherited from edm::stream::EDProducer<>
typedef CacheContexts< T...>	CacheTypes

typedef CacheTypes::GlobalCache	GlobalCache

typedef AbilityChecker< T...>	HasAbility

typedef CacheTypes::LuminosityBlockCache	LuminosityBlockCache

typedef LuminosityBlockContextT < LuminosityBlockCache, RunCache, GlobalCache >	LuminosityBlockContext

typedef CacheTypes::LuminosityBlockSummaryCache	LuminosityBlockSummaryCache

typedef CacheTypes::RunCache	RunCache

typedef RunContextT< RunCache, GlobalCache >	RunContext

typedef CacheTypes::RunSummaryCache	RunSummaryCache

Public Types inherited from edm::stream::EDProducerBase
typedef EDProducerAdaptorBase	ModuleType

Public Types inherited from edm::ProducerBase
typedef ProductRegistryHelper::TypeLabelList	TypeLabelList

Static Public Member Functions inherited from edm::stream::EDProducerBase
static const std::string &	baseType ()

static void	fillDescriptions (ConfigurationDescriptions &descriptions)

static void	prevalidate (ConfigurationDescriptions &descriptions)

Protected Member Functions inherited from edm::ProducerBase
void	callWhenNewProductsRegistered (std::function< void(BranchDescription const &)> const &func)

Protected Member Functions inherited from edm::EDConsumerBase
template<typename ProductType , BranchType B = InEvent>
EDGetTokenT< ProductType >	consumes (edm::InputTag const &tag)

EDGetToken	consumes (const TypeToGet &id, edm::InputTag const &tag)

template<BranchType B>
EDGetToken	consumes (TypeToGet const &id, edm::InputTag const &tag)

ConsumesCollector	consumesCollector ()
	Use a ConsumesCollector to gather consumes information from helper functions. More...

template<typename ProductType , BranchType B = InEvent>
void	consumesMany ()

void	consumesMany (const TypeToGet &id)

template<BranchType B>
void	consumesMany (const TypeToGet &id)

template<typename ProductType , BranchType B = InEvent>
EDGetTokenT< ProductType >	mayConsume (edm::InputTag const &tag)

EDGetToken	mayConsume (const TypeToGet &id, edm::InputTag const &tag)

template<BranchType B>
EDGetToken	mayConsume (const TypeToGet &id, edm::InputTag const &tag)

Detailed Description

Definition at line 29 of file TrackListMerger.cc.

Constructor & Destructor Documentation

TrackListMerger::TrackListMerger ( const edm::ParameterSet & conf )

explicit

Definition at line 199 of file TrackListMerger.cc.

References allowFirstHitShare_, copyExtras_, copyMVA_, edTokens(), epsilon_, edm::hlt::Exception, edm::ParameterSet::exists(), foundHitBonus_, edm::ParameterSet::getParameter(), edm::ParameterSet::getUntrackedParameter(), hasSelector_, i, indivShareFrac_, diffTwoXMLs::label, listsToMerge_, lostHitPenalty_, makeReKeyedSeeds_, maxNormalizedChisq_, minFound_, minPT_, convertSQLiteXML::ok, promoteQuality_, reco::TrackBase::qualityByName(), qualityToSet_, shareFrac_, AlCaHLTBitMon_QueryRunRegistry::string, trackProducers_, trkQualMod_, reco::TrackBase::undefQuality, and use_sharesInput_.

                                                             {
   copyExtras_ = conf.getUntrackedParameter<bool>("copyExtras", true);
   
   std::vector<edm::InputTag> trackProducerTags(conf.getParameter<std::vector<edm::InputTag> >("TrackProducers"));
   //which of these do I need to turn into vectors?
   maxNormalizedChisq_ =  conf.getParameter<double>("MaxNormalizedChisq");
   minPT_ =  conf.getParameter<double>("MinPT");
   minFound_ = (unsigned int)conf.getParameter<int>("MinFound");
   epsilon_ =  conf.getParameter<double>("Epsilon");
   shareFrac_ =  conf.getParameter<double>("ShareFrac");
   allowFirstHitShare_ = conf.getParameter<bool>("allowFirstHitShare");
   foundHitBonus_ = conf.getParameter<double>("FoundHitBonus");
   lostHitPenalty_ = conf.getParameter<double>("LostHitPenalty");
   indivShareFrac_=conf.getParameter<std::vector<double> >("indivShareFrac");
   std::string qualityStr = conf.getParameter<std::string>("newQuality");
   
   if (qualityStr != "") {
     qualityToSet_ = reco::TrackBase::qualityByName(conf.getParameter<std::string>("newQuality"));
   }
   else
     qualityToSet_ = reco::TrackBase::undefQuality;
   
   use_sharesInput_ = true;
   if ( epsilon_ > 0.0 )use_sharesInput_ = false;
   
   edm::VParameterSet setsToMerge=conf.getParameter<edm::VParameterSet>("setsToMerge");
   
   for ( unsigned int i=0; i<setsToMerge.size(); i++) {
     listsToMerge_.push_back(setsToMerge[i].getParameter<std::vector< int> >("tLists"));
     promoteQuality_.push_back(setsToMerge[i].getParameter<bool>("pQual"));
   }
   hasSelector_ = conf.getParameter<std::vector<int> >("hasSelector");
   copyMVA_     = conf.getParameter<bool>("copyMVA"); 
   
   std::vector<edm::InputTag> selectors(conf.getParameter<std::vector<edm::InputTag> >("selectedTrackQuals"));
   std::vector<edm::InputTag> mvaStores;
   if(conf.exists("mvaValueTags")){
     mvaStores = conf.getParameter<std::vector<edm::InputTag> >("mvaValueTags");
   }else{
     for (int i = 0; i < (int)selectors.size(); i++){
       edm::InputTag ntag(selectors[i].label(),"MVAVals");
       mvaStores.push_back(ntag);
     }
   }
   unsigned int numTrkColl = trackProducerTags.size();
   if (numTrkColl != hasSelector_.size() || numTrkColl != selectors.size()) {
     throw cms::Exception("Inconsistent size") << "need same number of track collections and selectors";
   }
   if (numTrkColl != hasSelector_.size() || numTrkColl != mvaStores.size()) {
     throw cms::Exception("Inconsistent size") << "need same number of track collections and MVA stores";
   }
   for (unsigned int i = indivShareFrac_.size(); i < numTrkColl; i++) {
     //      edm::LogWarning("TrackListMerger") << "No indivShareFrac for " << trackProducersTags <<". Using default value of 1";
     indivShareFrac_.push_back(1.0);
   }
   
   trkQualMod_=conf.getParameter<bool>("writeOnlyTrkQuals");
   if ( trkQualMod_) {
     bool ok=true;
     for ( unsigned int i=1; i<numTrkColl; i++) {
       if (!(trackProducerTags[i]==trackProducerTags[0])) ok=false;
     }
     if ( !ok) {
       throw cms::Exception("Bad input") << "to use writeOnlyTrkQuals=True all input InputTags must be the same";
     }
     produces<edm::ValueMap<int> >();
   }
   else{
     produces<reco::TrackCollection>();
     
     makeReKeyedSeeds_ = conf.getUntrackedParameter<bool>("makeReKeyedSeeds",false);
     if (makeReKeyedSeeds_){
       copyExtras_=true;
       produces<TrajectorySeedCollection>();
     }
     
     if (copyExtras_) {
       produces<reco::TrackExtraCollection>();
       produces<TrackingRecHitCollection>();
     }
     produces< std::vector<Trajectory> >();
     produces< TrajTrackAssociationCollection >();
   }
   produces<edm::ValueMap<float> >("MVAVals");
   
   // Do all the consumes
   trackProducers_.resize(numTrkColl);
   for (unsigned int i = 0; i < numTrkColl; ++i) {
     trackProducers_[i] = hasSelector_[i]>0 ? edTokens(trackProducerTags[i], selectors[i], mvaStores[i]) : edTokens(trackProducerTags[i], mvaStores[i]);
   }
 }

TrackListMerger::~TrackListMerger ( )

virtual

Definition at line 293 of file TrackListMerger.cc.

293 { }

Member Function Documentation

TkEDGetTokenss TrackListMerger::edTokens	(	const edm::InputTag &	tag,
		const edm::InputTag &	seltag,
		const edm::InputTag &	mvatag
	)

inlineprivate

Definition at line 69 of file TrackListMerger.cc.

Referenced by TrackListMerger().

                                                                                                           {
         return TkEDGetTokenss(tag, consumes<reco::TrackCollection>(tag), 
                                 consumes<std::vector<Trajectory> >(tag), consumes<TrajTrackAssociationCollection >(tag),
                                 consumes<edm::ValueMap<int> >(seltag), consumes<edm::ValueMap<float> >(mvatag));
     }

TkEDGetTokenss TrackListMerger::edTokens	(	const edm::InputTag &	tag,
		const edm::InputTag &	mvatag
	)

inlineprivate

Definition at line 74 of file TrackListMerger.cc.

                                                                                {
         return TkEDGetTokenss(tag, consumes<reco::TrackCollection>(tag), 
                                 consumes<std::vector<Trajectory> >(tag), consumes<TrajTrackAssociationCollection >(tag),
                                 edm::EDGetTokenT<edm::ValueMap<int> >(), consumes<edm::ValueMap<float> >(mvatag));
     }

void TrackListMerger::produce	(	edm::Event &	e,
		const edm::EventSetup &	c
	)

overridevirtual

Implements edm::stream::EDProducerBase.

Definition at line 296 of file TrackListMerger.cc.

References funct::abs(), reco::TrackExtraBase::add(), reco::TrackBase::algo(), allowFirstHitShare_, asciidump::at, edm::OwnVector< T, P >::back(), begin, reco::TrackBase::chi2(), edm::AssociationMap< edm::OneToOne< std::vector< Trajectory >, reco::TrackCollection, unsigned short > >::const_iterator, constexpr, copyExtras_, copyMVA_, delta, TrajectorySeed::direction(), end, epsilon_, edm::HandleBase::failedToGet(), spr::find(), first, foundHitBonus_, TrackingRecHit::geographicalId(), edm::Event::getByLabel(), edm::Event::getByToken(), edm::Event::getProvenance(), edm::Event::getRefBeforePut(), hasSelector_, reco::TrackBase::highPurity, reco::TrackBase::hitPattern(), i, edm::RefProd< T >::id(), edm::Ref< C, T, F >::id(), cuy::ii, indivShareFrac_, reco::Track::innerDetId(), reco::Track::innerMomentum(), reco::Track::innerOk(), reco::Track::innerPosition(), reco::Track::innerStateCovariance(), TrackingRecHit::isValid(), j, findQualityFiles::jj, edm::Ref< C, T, F >::key(), likely, listsToMerge_, TrackingRecHit::localPosition(), reco::TrackBase::loose, lostHitPenalty_, makeReKeyedSeeds_, match(), maxNormalizedChisq_, bookConverter::min, minFound_, minPT_, edm::Provenance::moduleLabel(), n, TrajectorySeed::nHits(), reco::TrackBase::numberOfLostHits(), reco::TrackBase::numberOfValidHits(), reco::HitPattern::numberOfValidPixelHits(), reco::Track::outerDetId(), reco::Track::outerMomentum(), reco::Track::outerOk(), reco::Track::outerPosition(), reco::Track::outerStateCovariance(), outputSeeds, outputTrajs, outputTrkExtras, outputTrkHits, outputTrks, outputTTAss, edm::Provenance::processName(), edm::Provenance::productInstanceName(), promoteQuality_, edm::OwnVector< T, P >::push_back(), edm::Event::put(), qualityToSet_, TrackingRecHit::rawId(), reco::Track::recHit(), TrajectorySeed::recHits(), reco::Track::recHitsBegin(), reco::Track::recHitsEnd(), reco::Track::recHitsSize(), refTrajs, refTrajSeeds, refTrkExtras, refTrkHits, refTrks, ClusterRemovalRefSetter::reKey(), edm::OwnVector< T, P >::reserve(), reco::Track::residuals(), reco::Track::seedDirection(), reco::Track::seedRef(), reco::TrackExtra::setResiduals(), shareFrac_, TrackingRecHit::some, TrajectorySeed::startingState(), DetId::subdetId(), reco::TrackBase::tight, trackProducers_, trkQualMod_, unlikely, use_sharesInput_, and PV3DBase< T, PVType, FrameType >::x().

 {
   // extract tracker geometry
   //
   //edm::ESHandle<TrackerGeometry> theG;
   //es.get<TrackerDigiGeometryRecord>().get(theG);
   
   //    using namespace reco;
   
   // get Inputs
   // if 1 input list doesn't exist, make an empty list, issue a warning, and continue
   // this allows TrackListMerger to be used as a cleaner only if handed just one list
   // if both input lists don't exist, will issue 2 warnings and generate an empty output collection
   //
   static const reco::TrackCollection s_empty;
   
   std::vector<const reco::TrackCollection *> trackColls;
   std::vector<edm::Handle<reco::TrackCollection> > trackHandles(trackProducers_.size());
   for ( unsigned int i=0; i<trackProducers_.size(); i++) {
     trackColls.push_back(0);
     //edm::Handle<reco::TrackCollection> trackColl;
     e.getByToken(trackProducers_[i].tk, trackHandles[i]);
     if (trackHandles[i].isValid()) {
       trackColls[i]= trackHandles[i].product();
     } else {
       edm::LogWarning("TrackListMerger") << "TrackCollection " << trackProducers_[i].tag <<" not found";
       trackColls[i]=&s_empty;
     }
   }
   
   unsigned int collsSize =trackColls.size();
   unsigned int rSize=0;
   unsigned int trackCollSizes[collsSize];
   unsigned int trackCollFirsts[collsSize];
   for (unsigned int i=0; i!=collsSize; i++) {
     trackCollSizes[i]=trackColls[i]->size();
     trackCollFirsts[i]=rSize;
     rSize+=trackCollSizes[i];
   }
   
   statCount.begin(rSize);
   
   //
   //  quality cuts first
   //
   int i=-1;
   
   int selected[rSize];
   int indexG[rSize];
   bool trkUpdated[rSize];
   int trackCollNum[rSize];
   int trackQuals[rSize];
   float trackMVAs[rSize];
   for (unsigned int j=0; j<rSize;j++) {
     indexG[j]=-1; selected[j]=1; trkUpdated[j]=false; trackCollNum[j]=0; trackQuals[j]=0;trackMVAs[j] = -99.0;
   }
   
   int ngood=0;
   for (unsigned int j=0; j!= collsSize; j++) {
     const reco::TrackCollection *tC1=trackColls[j];
     
     edm::Handle<edm::ValueMap<int> > trackSelColl;
     edm::Handle<edm::ValueMap<float> > trackMVAStore;
     if ( copyMVA_ )
       e.getByToken(trackProducers_[j].tmva, trackMVAStore);
     if ( hasSelector_[j]>0 ){
       e.getByToken(trackProducers_[j].tsel, trackSelColl);
     }
     
     if ( 0<tC1->size() ){
       unsigned int iC=0;
       for (reco::TrackCollection::const_iterator track=tC1->begin(); track!=tC1->end(); track++){
     i++;
     trackCollNum[i]=j;
     trackQuals[i]=track->qualityMask();
     
     reco::TrackRef trkRef=reco::TrackRef(trackHandles[j],iC);
     if ( copyMVA_ )
       if( (*trackMVAStore).contains(trkRef.id()) ) trackMVAs[i] = (*trackMVAStore)[trkRef];
     if ( hasSelector_[j]>0 ) {
       int qual=(*trackSelColl)[trkRef];
       if ( qual < 0 ) {
         selected[i]=0;
         iC++;
         continue;
       }else{
         trackQuals[i]=qual;
       }
     }
     iC++;
     selected[i]=trackQuals[i]+10;//10 is magic number used throughout...
     if ((short unsigned)track->ndof() < 1){
       selected[i]=0;
       continue;
     }
     if (track->normalizedChi2() > maxNormalizedChisq_){
       selected[i]=0;
       continue;
     }
     if (track->found() < minFound_){
       selected[i]=0;
       continue;
     }
     if (track->pt() < minPT_){
       selected[i]=0;
       continue;
     }
     // good!
     indexG[i] = ngood++;
     //if ( beVerb) std::cout << "inverb " << track->pt() << " " << selected[i] << std::endl;
       }//end loop over tracks
     }//end more than 0 track
   } // loop over trackcolls
   
   
   statCount.pre(ngood);
   
   //cache the id and rechits of valid hits
   typedef std::pair<unsigned int, const TrackingRecHit*> IHit;
   std::vector<IHit> rh1[ngood];  // an array of vectors!
   //const TrackingRecHit*  fh1[ngood];  // first hit...
   uint8_t algo[ngood];
   float score[ngood];
   
   
   for ( unsigned int j=0; j<rSize; j++) {
     if (selected[j]==0) continue;
     int i = indexG[j];
     assert(i>=0);
     unsigned int collNum=trackCollNum[j];
     unsigned int trackNum=j-trackCollFirsts[collNum];
     const reco::Track *track=&((trackColls[collNum])->at(trackNum));
     
     algo[i]=track->algo();
     int validHits=track->numberOfValidHits();
     int validPixelHits=track->hitPattern().numberOfValidPixelHits();
     int lostHits=track->numberOfLostHits();
     score[i] = foundHitBonus_*validPixelHits+foundHitBonus_*validHits - lostHitPenalty_*lostHits - track->chi2();
     
     
     rh1[i].reserve(validHits) ;
     auto compById = [](IHit const &  h1, IHit const & h2) {return h1.first < h2.first;};
     // fh1[i] = &(**track->recHitsBegin());
     for (trackingRecHit_iterator it = track->recHitsBegin();  it != track->recHitsEnd(); ++it) {
       const TrackingRecHit* hit = &(**it);
       unsigned int id = hit->rawId() ;
       if(hit->geographicalId().subdetId()>2)  id &= (~3); // mask mono/stereo in strips...
       if likely(hit->isValid()) { rh1[i].emplace_back(id,hit); std::push_heap(rh1[i].begin(),rh1[i].end(),compById); }
     }
     std::sort_heap(rh1[i].begin(),rh1[i].end(),compById);
   }
   
   //DL here
   if likely(ngood>1 && collsSize>1)
          for ( unsigned int ltm=0; ltm<listsToMerge_.size(); ltm++) {
            int saveSelected[rSize];
            bool notActive[collsSize];
            for (unsigned int cn=0;cn!=collsSize;++cn)
          notActive[cn]= find(listsToMerge_[ltm].begin(),listsToMerge_[ltm].end(),cn)==listsToMerge_[ltm].end();
            
            for ( unsigned int i=0; i<rSize; i++) saveSelected[i]=selected[i];
            
            //DL protect against 0 tracks?
            for ( unsigned int i=0; i<rSize-1; i++) {
          if (selected[i]==0) continue;
          unsigned int collNum=trackCollNum[i];
          
          //check that this track is in one of the lists for this iteration
          if (notActive[collNum]) continue;
          
          int k1 = indexG[i];
          unsigned int nh1=rh1[k1].size();
          int qualityMaskT1 = trackQuals[i];
          
          int nhit1 = nh1; // validHits[k1];
          float score1 = score[k1];
          
          // start at next collection
          for ( unsigned int j=i+1; j<rSize; j++) {
            if (selected[j]==0) continue;
            unsigned int collNum2=trackCollNum[j];
            if ( (collNum == collNum2) && indivShareFrac_[collNum] > 0.99) continue;
            //check that this track is in one of the lists for this iteration
            if (notActive[collNum2]) continue;
            
            int k2 = indexG[j];
            
            int newQualityMask = -9; //avoid resetting quality mask if not desired 10+ -9 =1
            if (promoteQuality_[ltm]) {
              int maskT1= saveSelected[i]>1? saveSelected[i]-10 : qualityMaskT1;
              int maskT2= saveSelected[j]>1? saveSelected[j]-10 : trackQuals[j];
              newQualityMask =(maskT1 | maskT2); // take OR of trackQuality
            }
            unsigned int nh2=rh1[k2].size();
            int nhit2 = nh2;
            
            
            auto share = use_sharesInput_ ?
              [](const TrackingRecHit*  it,const TrackingRecHit*  jt, float)->bool { return it->sharesInput(jt,TrackingRecHit::some); } :
            [](const TrackingRecHit*  it,const TrackingRecHit*  jt, float eps)->bool {
              float delta = std::abs ( it->localPosition().x()-jt->localPosition().x() );
              return (it->geographicalId()==jt->geographicalId())&&(delta<eps);
            };
            
            statCount.start();
            
            //loop over rechits
            int noverlap=0;
            int firstoverlap=0;
            // check first hit  (should use REAL first hit?)
            if unlikely(allowFirstHitShare_ && rh1[k1][0].first==rh1[k2][0].first ) {
                const TrackingRecHit*  it = rh1[k1][0].second;
                const TrackingRecHit*  jt = rh1[k2][0].second;
                if (share(it,jt,epsilon_)) firstoverlap=1;
              }
            
            
            // exploit sorting
            unsigned int jh=0;
            unsigned int ih=0;
            while (ih!=nh1 && jh!=nh2) {
              // break if not enough to go...
              // if ( nprecut-noverlap+firstoverlap > int(nh1-ih)) break;
              // if ( nprecut-noverlap+firstoverlap > int(nh2-jh)) break;
              auto const id1 = rh1[k1][ih].first;
              auto const id2 = rh1[k2][jh].first;
              if (id1<id2) ++ih;
              else if (id2<id1) ++jh;
              else {
                // in case of split-hit do full conbinatorics
                auto li=ih; while( (++li)!=nh1 && id1 == rh1[k1][li].first);
                auto lj=jh; while( (++lj)!=nh2 && id2 == rh1[k2][lj].first);
                for (auto ii=ih; ii!=li; ++ii)
              for (auto jj=jh; jj!=lj; ++jj) {
                const TrackingRecHit*  it = rh1[k1][ii].second;
                const TrackingRecHit*  jt = rh1[k2][jj].second;
                if (share(it,jt,epsilon_))  noverlap++;
              }
                jh=lj; ih=li;
              } // equal ids
              
            } //loop over ih & jh
            
            bool dupfound = (collNum != collNum2) ? (noverlap-firstoverlap) > (std::min(nhit1,nhit2)-firstoverlap)*shareFrac_ :
              (noverlap-firstoverlap) > (std::min(nhit1,nhit2)-firstoverlap)*indivShareFrac_[collNum];
            
            if ( dupfound ) {
              float score2 = score[k2];
              constexpr float almostSame = 0.01f; // difference rather than ratio due to possible negative values for score
              if ( score1 - score2 > almostSame ) {
                selected[j]=0;
                selected[i]=10+newQualityMask; // add 10 to avoid the case where mask = 1
                trkUpdated[i]=true;
              } else if ( score2 - score1 > almostSame ) {
                selected[i]=0;
                selected[j]=10+newQualityMask;  // add 10 to avoid the case where mask = 1
                trkUpdated[j]=true;
              }else{
                // If tracks from both iterations are virtually identical, choose the one with the best quality or with lower algo
                if ((trackQuals[j] & (1<<reco::TrackBase::loose|1<<reco::TrackBase::tight|1<<reco::TrackBase::highPurity) ) ==
                (trackQuals[i] & (1<<reco::TrackBase::loose|1<<reco::TrackBase::tight|1<<reco::TrackBase::highPurity) )) {
              //same quality, pick earlier algo
              if (algo[k1] <= algo[k2]) {
                selected[j]=0;
                selected[i]=10+newQualityMask; // add 10 to avoid the case where mask = 1
                trkUpdated[i]=true;
              } else {
                selected[i]=0;
                selected[j]=10+newQualityMask; // add 10 to avoid the case where mask = 1
                trkUpdated[j]=true;        
              }
                } else if ((trackQuals[j] & (1<<reco::TrackBase::loose|1<<reco::TrackBase::tight|1<<reco::TrackBase::highPurity) ) <
                   (trackQuals[i] & (1<<reco::TrackBase::loose|1<<reco::TrackBase::tight|1<<reco::TrackBase::highPurity) )) {
              selected[j]=0;
              selected[i]=10+newQualityMask; // add 10 to avoid the case where mask = 1
              trkUpdated[i]=true;
                }else{
              selected[i]=0;
              selected[j]=10+newQualityMask; // add 10 to avoid the case where mask = 1
              trkUpdated[j]=true;
                }
              }//end fi < fj
              statCount.overlap();
              /*
                if (at0[k1]&&at0[k2]) {
                statCount.dp(dphi);
                if (dz<1.f) statCount.de(deta);
                }
              */
            }//end got a duplicate
            else {
              statCount.noOverlap();
            }
            //stop if the ith track is now unselected
            if (selected[i]==0) break;
          }//end track2 loop
            }//end track loop
          } //end loop over track list sets
   
   
   
   std::auto_ptr<edm::ValueMap<float> > vmMVA = std::auto_ptr<edm::ValueMap<float> >(new edm::ValueMap<float>);
   edm::ValueMap<float>::Filler fillerMVA(*vmMVA);
   
   
   
   // special case - if just doing the trkquals
   if (trkQualMod_) {
     unsigned int tSize=trackColls[0]->size();
     std::auto_ptr<edm::ValueMap<int> > vm = std::auto_ptr<edm::ValueMap<int> >(new edm::ValueMap<int>);
     edm::ValueMap<int>::Filler filler(*vm);
     
     std::vector<int> finalQuals(tSize,-1); //default is unselected
     for ( unsigned int i=0; i<rSize; i++) {
       unsigned int tNum=i%tSize;
       
       if (selected[i]>1 ) {
     finalQuals[tNum]=selected[i]-10;
     if (trkUpdated[i])
       finalQuals[tNum]=(finalQuals[tNum] | (1<<qualityToSet_));
       }
       if ( selected[i]==1 )
     finalQuals[tNum]=trackQuals[i];
     }
     
     filler.insert(trackHandles[0], finalQuals.begin(),finalQuals.end());
     filler.fill();
     
     e.put(vm);
     
     std::vector<float> mvaVec(tSize,-99);
     
     for(unsigned int i = 0; i < rSize; i++){
       unsigned int tNum = i % tSize;
       mvaVec[tNum] = trackMVAs[tNum];
     }
     
     fillerMVA.insert(trackHandles[0],mvaVec.begin(),mvaVec.end());
     fillerMVA.fill();
     if ( copyMVA_) 
       e.put(vmMVA,"MVAVals");
     return;
   }
   
   
   //
   //  output selected tracks - if any
   //
   
   reco::TrackRef trackRefs[rSize];
   edm::RefToBase<TrajectorySeed> seedsRefs[rSize];
   
   unsigned int nToWrite=0;
   for ( unsigned int i=0; i<rSize; i++)
     if (selected[i]!=0) nToWrite++;
   
   std::vector<float> mvaVec;
   
   outputTrks = std::auto_ptr<reco::TrackCollection>(new reco::TrackCollection);
   outputTrks->reserve(nToWrite);
   refTrks = e.getRefBeforePut<reco::TrackCollection>();
   
   if (copyExtras_) {
     outputTrkExtras = std::auto_ptr<reco::TrackExtraCollection>(new reco::TrackExtraCollection);
     outputTrkExtras->reserve(nToWrite);
     refTrkExtras    = e.getRefBeforePut<reco::TrackExtraCollection>();
     outputTrkHits   = std::auto_ptr<TrackingRecHitCollection>(new TrackingRecHitCollection);
     outputTrkHits->reserve(nToWrite*25);
     refTrkHits      = e.getRefBeforePut<TrackingRecHitCollection>();
     if (makeReKeyedSeeds_){
       outputSeeds = std::auto_ptr<TrajectorySeedCollection>(new TrajectorySeedCollection);
       outputSeeds->reserve(nToWrite);
       refTrajSeeds = e.getRefBeforePut<TrajectorySeedCollection>();
     }
   }
   
   
   outputTrajs = std::auto_ptr< std::vector<Trajectory> >(new std::vector<Trajectory>());
   outputTrajs->reserve(rSize);
   outputTTAss = std::auto_ptr< TrajTrackAssociationCollection >(new TrajTrackAssociationCollection());
   
   
   
   for ( unsigned int i=0; i<rSize; i++) {
     if (selected[i]==0) {
       trackRefs[i]=reco::TrackRef();
       continue;
     }
     
     unsigned int collNum=trackCollNum[i];
     unsigned int trackNum=i-trackCollFirsts[collNum];
     const reco::Track *track=&((trackColls[collNum])->at(trackNum));
     outputTrks->push_back( reco::Track( *track ) );
     mvaVec.push_back(trackMVAs[i]);
     if (selected[i]>1 ) {
       outputTrks->back().setQualityMask(selected[i]-10);
       if (trkUpdated[i]) 
     outputTrks->back().setQuality(qualityToSet_);
     }
     //might duplicate things, but doesnt hurt
     if ( selected[i]==1 ) 
       outputTrks->back().setQualityMask(trackQuals[i]);
     
     // if ( beVerb ) std::cout << "selected " << outputTrks->back().pt() << " " << outputTrks->back().qualityMask() << " " << selected[i] << std::endl;
     
     //fill the TrackCollection
     if (copyExtras_) {
       edm::RefToBase<TrajectorySeed> origSeedRef = track->seedRef();
       //creating a seed with rekeyed clusters if required
       if (makeReKeyedSeeds_){
     bool doRekeyOnThisSeed=false;
     
     edm::InputTag clusterRemovalInfos("");
     //grab on of the hits of the seed
     if (origSeedRef->nHits()!=0){
       TrajectorySeed::const_iterator firstHit=origSeedRef->recHits().first;
       const TrackingRecHit *hit = &*firstHit;
       if (firstHit->isValid()){
         edm::ProductID  pID=clusterProductB(hit);
         // the cluster collection either produced a removalInfo or mot
         //get the clusterremoval info from the provenance: will rekey if this is found
         edm::Handle<reco::ClusterRemovalInfo> CRIh;
         edm::Provenance prov=e.getProvenance(pID);
         clusterRemovalInfos=edm::InputTag(prov.moduleLabel(),
                           prov.productInstanceName(),
                           prov.processName());
         doRekeyOnThisSeed=e.getByLabel(clusterRemovalInfos,CRIh);
       }//valid hit
     }//nhit!=0
     
     if (doRekeyOnThisSeed && !(clusterRemovalInfos==edm::InputTag(""))) {
       ClusterRemovalRefSetter refSetter(e,clusterRemovalInfos);
       TrajectorySeed::recHitContainer  newRecHitContainer;
       newRecHitContainer.reserve(origSeedRef->nHits());
       TrajectorySeed::const_iterator iH=origSeedRef->recHits().first;
       TrajectorySeed::const_iterator iH_end=origSeedRef->recHits().second;
       for (;iH!=iH_end;++iH){
         newRecHitContainer.push_back(*iH);
         refSetter.reKey(&newRecHitContainer.back());
       }
       outputSeeds->push_back( TrajectorySeed( origSeedRef->startingState(),
                           newRecHitContainer,
                           origSeedRef->direction()));
     }
     //doRekeyOnThisSeed=true
     else{
       //just copy the one we had before
       outputSeeds->push_back( TrajectorySeed(*origSeedRef));
     }
     edm::Ref<TrajectorySeedCollection> pureRef(refTrajSeeds, outputSeeds->size()-1);
     origSeedRef=edm::RefToBase<TrajectorySeed>( pureRef);
       }//creating a new seed and rekeying it rechit clusters.
       
       // Fill TrackExtra collection
       outputTrkExtras->push_back( reco::TrackExtra(
                            track->outerPosition(), track->outerMomentum(), track->outerOk(),
                            track->innerPosition(), track->innerMomentum(), track->innerOk(),
                            track->outerStateCovariance(), track->outerDetId(),
                            track->innerStateCovariance(), track->innerDetId(),
                            track->seedDirection(), origSeedRef ) );
       seedsRefs[i]=origSeedRef;
       outputTrks->back().setExtra( reco::TrackExtraRef( refTrkExtras, outputTrkExtras->size() - 1) );
       reco::TrackExtra & tx = outputTrkExtras->back();
       tx.setResiduals(track->residuals());
       
       // fill TrackingRecHits
       unsigned nh1=track->recHitsSize();
       for ( unsigned ih=0; ih<nh1; ++ih ) {
     //const TrackingRecHit*hit=&((*(track->recHit(ih))));
     outputTrkHits->push_back( track->recHit(ih)->clone() );
     tx.add( TrackingRecHitRef( refTrkHits, outputTrkHits->size() - 1) );
       }
     }
     trackRefs[i] = reco::TrackRef(refTrks, outputTrks->size() - 1);
     
     
   }//end faux loop over tracks
   
   //Fill the trajectories, etc. for 1st collection
   refTrajs    = e.getRefBeforePut< std::vector<Trajectory> >();
   
   for (unsigned int ti=0; ti<trackColls.size(); ti++) {
     edm::Handle< std::vector<Trajectory> >  hTraj1;
     edm::Handle< TrajTrackAssociationCollection >  hTTAss1;
     e.getByToken(trackProducers_[ti].traj, hTraj1);
     e.getByToken(trackProducers_[ti].tass, hTTAss1);
     
     if (hTraj1.failedToGet() || hTTAss1.failedToGet()) continue;
     
     for (size_t i = 0, n = hTraj1->size(); i < n; ++i) {
       edm::Ref< std::vector<Trajectory> > trajRef(hTraj1, i);
       TrajTrackAssociationCollection::const_iterator match = hTTAss1->find(trajRef);
       if (match != hTTAss1->end()) {
     const edm::Ref<reco::TrackCollection> &trkRef = match->val;
     short oldKey = trackCollFirsts[ti]+static_cast<short>(trkRef.key());
     if (trackRefs[oldKey].isNonnull()) {
       outputTrajs->push_back( *trajRef );
       //if making extras and the seeds at the same time, change the seed ref on the trajectory
       if (copyExtras_ && makeReKeyedSeeds_)
         outputTrajs->back().setSeedRef( seedsRefs[oldKey] );
       outputTTAss->insert ( edm::Ref< std::vector<Trajectory> >(refTrajs, outputTrajs->size() - 1),
                 trackRefs[oldKey] );
     }
       }
     }
   }
   
   statCount.end(outputTrks->size());
   
   edm::ProductID nPID = refTrks.id();
   edm::TestHandle<reco::TrackCollection> outHandle(outputTrks.get(),nPID);
   fillerMVA.insert(outHandle,mvaVec.begin(),mvaVec.end());
   fillerMVA.fill();
   
   e.put(outputTrks);
   if ( copyMVA_ )
     e.put(vmMVA,"MVAVals");
   if (copyExtras_) {
     e.put(outputTrkExtras);
     e.put(outputTrkHits);
     if (makeReKeyedSeeds_)
       e.put(outputSeeds);
   }
   e.put(outputTrajs);
   e.put(outputTTAss);
   return;
   
 }//end produce