/data/doxygen/doxygen-1.7.3/gen/CMSSW_4_2_8/src/RecoTracker/FinalTrackSelectors/src/SimpleTrackListMerger.cc

Go to the documentation of this file.

//
// Package:         RecoTracker/FinalTrackSelectors
// Class:           SimpleTrackListMerger
// 
// Description:     TrackList Cleaner and Merger
//
// Original Author: Steve Wagner, stevew@pizero.colorado.edu
// Created:         Sat Jan 14 22:00:00 UTC 2006
//
// $Author: innocent $
// $Date: 2010/12/14 15:03:08 $
// $Revision: 1.27 $
//

#include <memory>
#include <string>
#include <iostream>
#include <cmath>
#include <vector>

#include "RecoTracker/FinalTrackSelectors/interface/SimpleTrackListMerger.h"

#include "DataFormats/TrackerRecHit2D/interface/SiStripMatchedRecHit2DCollection.h"
#include "DataFormats/TrackerRecHit2D/interface/SiStripRecHit2DCollection.h"
#include "DataFormats/TrackerRecHit2D/interface/SiStripRecHit1DCollection.h"
#include "DataFormats/TrajectorySeed/interface/TrajectorySeedCollection.h"
#include "DataFormats/TrackCandidate/interface/TrackCandidateCollection.h"

#include "FWCore/Framework/interface/ESHandle.h"

#include "FWCore/MessageLogger/interface/MessageLogger.h"

#include "Geometry/CommonDetUnit/interface/GeomDetUnit.h"
#include "Geometry/CommonDetUnit/interface/GeomDet.h"
#include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
#include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"

//#include "DataFormats/TrackReco/src/classes.h"

#include "RecoTracker/TrackProducer/interface/ClusterRemovalRefSetter.h"


namespace cms
{
  
  edm::ProductID clusterProduct( const TrackingRecHit *hit){
    edm::ProductID pID;
    //cast it into the proper class     and find productID
    DetId detid = hit->geographicalId(); 
    uint32_t subdet = detid.subdetId();
    if ((subdet == PixelSubdetector::PixelBarrel) || (subdet == PixelSubdetector::PixelEndcap)) {
      pID=reinterpret_cast<const SiPixelRecHit *>(hit)->cluster().id();
    } else {
      const std::type_info &type = typeid(*hit);
      if (type == typeid(SiStripRecHit2D)) {
        pID=reinterpret_cast<const SiStripRecHit2D *>(hit)->cluster().id();
      } else if (type == typeid(SiStripRecHit1D)) {
        pID=reinterpret_cast<const SiStripRecHit1D *>(hit)->cluster().id();     
      } else if (type == typeid(SiStripMatchedRecHit2D)) {
        const SiStripMatchedRecHit2D *mhit = reinterpret_cast<const SiStripMatchedRecHit2D *>(hit);
        pID=mhit->monoHit()->cluster().id();
      } else if (type == typeid(ProjectedSiStripRecHit2D)) {
        const ProjectedSiStripRecHit2D *phit = reinterpret_cast<const ProjectedSiStripRecHit2D *>(hit);
        pID=(&phit->originalHit())->cluster().id();
      } else throw cms::Exception("Unknown RecHit Type") << "RecHit of type " << type.name() << " not supported. (use c++filt to demangle the name)";
    }
        
    return pID;}
  


  SimpleTrackListMerger::SimpleTrackListMerger(edm::ParameterSet const& conf) : 
    conf_(conf)
  {
    copyExtras_ = conf_.getUntrackedParameter<bool>("copyExtras", true);

    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 >();
    

  }


  // Virtual destructor needed.
  SimpleTrackListMerger::~SimpleTrackListMerger() { }  

  // Functions that gets called by framework every event
  void SimpleTrackListMerger::produce(edm::Event& e, const edm::EventSetup& es)
  {
    // retrieve producer name of input TrackCollection(s)
    std::string trackProducer1 = conf_.getParameter<std::string>("TrackProducer1");
    std::string trackProducer2 = conf_.getParameter<std::string>("TrackProducer2");

    double maxNormalizedChisq =  conf_.getParameter<double>("MaxNormalizedChisq");
    double minPT =  conf_.getParameter<double>("MinPT");
    unsigned int minFound = (unsigned int)conf_.getParameter<int>("MinFound");
    double epsilon =  conf_.getParameter<double>("Epsilon");
    bool use_sharesInput = true;
    if ( epsilon > 0.0 )use_sharesInput = false;
    double shareFrac =  conf_.getParameter<double>("ShareFrac");
  
    bool promoteQuality = conf_.getParameter<bool>("promoteTrackQuality");
    bool allowFirstHitShare = conf_.getParameter<bool>("allowFirstHitShare");
//

    // New track quality should be read from the file
    std::string qualityStr = conf_.getParameter<std::string>("newQuality");
    reco::TrackBase::TrackQuality qualityToSet;
    if (qualityStr != "") {
      qualityToSet = reco::TrackBase::qualityByName(conf_.getParameter<std::string>("newQuality"));
    }
    else 
      qualityToSet = reco::TrackBase::undefQuality;

    // 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 SimpleTrackListMerger 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
    //
    const reco::TrackCollection *TC1 = 0;
    static const reco::TrackCollection s_empty1, s_empty2;
    edm::Handle<reco::TrackCollection> trackCollection1;
    e.getByLabel(trackProducer1, trackCollection1);
    if (trackCollection1.isValid()) {
      TC1 = trackCollection1.product();
      //std::cout << "1st collection " << trackProducer1 << " has "<< TC1->size() << " tracks" << std::endl ;
    } else {
      TC1 = &s_empty1;
      edm::LogWarning("SimpleTrackListMerger") << "1st TrackCollection " << trackProducer1 << " not found; will only clean 2nd TrackCollection " << trackProducer2 ;
    }
    const reco::TrackCollection tC1 = *TC1;

    const reco::TrackCollection *TC2 = 0;
    edm::Handle<reco::TrackCollection> trackCollection2;
    e.getByLabel(trackProducer2, trackCollection2);
    if (trackCollection2.isValid()) {
      TC2 = trackCollection2.product();
      //std::cout << "2nd collection " << trackProducer2 << " has "<< TC2->size() << " tracks" << std::endl ;
    } else {
        TC2 = &s_empty2;
        edm::LogWarning("SimpleTrackListMerger") << "2nd TrackCollection " << trackProducer2 << " not found; will only clean 1st TrackCollection " << trackProducer1 ;
    }
    const reco::TrackCollection tC2 = *TC2;

    // Step B: create empty output collection
    outputTrks = std::auto_ptr<reco::TrackCollection>(new reco::TrackCollection);
    refTrks = e.getRefBeforePut<reco::TrackCollection>();      

    if (copyExtras_) {
        outputTrkExtras = std::auto_ptr<reco::TrackExtraCollection>(new reco::TrackExtraCollection);
        outputTrkExtras->reserve(TC1->size()+TC2->size());
        refTrkExtras    = e.getRefBeforePut<reco::TrackExtraCollection>();
        outputTrkHits   = std::auto_ptr<TrackingRecHitCollection>(new TrackingRecHitCollection);
        outputTrkHits->reserve((TC1->size()+TC2->size())*25);
        refTrkHits      = e.getRefBeforePut<TrackingRecHitCollection>();
        if (makeReKeyedSeeds_){
          outputSeeds = std::auto_ptr<TrajectorySeedCollection>(new TrajectorySeedCollection);
          outputSeeds->reserve(TC1->size()+TC2->size());
          refTrajSeeds = e.getRefBeforePut<TrajectorySeedCollection>();
        }
    }

    outputTrajs = std::auto_ptr< std::vector<Trajectory> >(new std::vector<Trajectory>()); 
    outputTrajs->reserve(TC1->size()+TC2->size());
    outputTTAss = std::auto_ptr< TrajTrackAssociationCollection >(new TrajTrackAssociationCollection());
    //outputTTAss->reserve(TC1->size()+TC2->size());//how do I reserve space for an association map?

  //
  //  no input tracks
  //

//    if ( tC1.empty() ){
//      LogDebug("RoadSearch") << "Found " << output.size() << " clouds.";
//      e.put(output);
//      return;  
//    }

  //
  //  quality cuts first
  // 
    int i;

    std::vector<int> selected1; for (unsigned int i=0; i<tC1.size(); ++i){selected1.push_back(1);}

   if ( 0<tC1.size() ){
      i=-1;
      for (reco::TrackCollection::const_iterator track=tC1.begin(); track!=tC1.end(); track++){
        i++;
        if ((short unsigned)track->ndof() < 1){
          selected1[i]=0; 
          //std::cout << "L1Track "<< i << " rejected in SimpleTrackListMerger; ndof() < 1" << std::endl ;
          continue;
        }
        if (track->normalizedChi2() > maxNormalizedChisq){
          selected1[i]=0; 
          //std::cout << "L1Track "<< i << " rejected in SimpleTrackListMerger; normalizedChi2() > maxNormalizedChisq " << track->normalizedChi2() << " " << maxNormalizedChisq << std::endl ;
          continue;
        }
        if (track->found() < minFound){
          selected1[i]=0; 
          //std::cout << "L1Track "<< i << " rejected in SimpleTrackListMerger; found() < minFound " << track->found() << " " << minFound << std::endl ;
          continue;
        }
        if (track->pt() < minPT){
          selected1[i]=0; 
          //std::cout << "L1Track "<< i << " rejected in SimpleTrackListMerger; pt() < minPT " << track->pt() << " " << minPT << std::endl ;
          continue;
        }
      }//end loop over tracks
   }//end more than 0 track


    std::vector<int> selected2; for (unsigned int i=0; i<tC2.size(); ++i){selected2.push_back(1);}

   if ( 0<tC2.size() ){
      i=-1;
      for (reco::TrackCollection::const_iterator track=tC2.begin(); track!=tC2.end(); track++){
        i++;
        if ((short unsigned)track->ndof() < 1){
          selected2[i]=0; 
          //std::cout << "L2Track "<< i << " rejected in SimpleTrackListMerger; ndof() < 1" << std::endl ;
          continue;
        }
        if (track->normalizedChi2() > maxNormalizedChisq){
          selected2[i]=0; 
          //std::cout << "L2Track "<< i << " rejected in SimpleTrackListMerger; normalizedChi2() > maxNormalizedChisq " << track->normalizedChi2() << " " << maxNormalizedChisq << std::endl ;
          continue;
        }
        if (track->found() < minFound){
          selected2[i]=0; 
          //std::cout << "L2Track "<< i << " rejected in SimpleTrackListMerger; found() < minFound " << track->found() << " " << minFound << std::endl ;
          continue;
        }
        if (track->pt() < minPT){
          selected2[i]=0; 
          //std::cout << "L2Track "<< i << " rejected in SimpleTrackListMerger; pt() < minPT " << track->pt() << " " << minPT << std::endl ;
          continue;
        }
      }//end loop over tracks
   }//end more than 0 track

   std::map<reco::TrackCollection::const_iterator, std::vector<const TrackingRecHit*> > rh1;
   std::map<reco::TrackCollection::const_iterator, std::vector<const TrackingRecHit*> > rh2;
   for (reco::TrackCollection::const_iterator track=tC1.begin(); track!=tC1.end(); ++track){
     trackingRecHit_iterator itB = track->recHitsBegin();
     trackingRecHit_iterator itE = track->recHitsEnd();
     for (trackingRecHit_iterator it = itB;  it != itE; ++it) { 
       const TrackingRecHit* hit = &(**it);
       rh1[track].push_back(hit);
     }
   }
   for (reco::TrackCollection::const_iterator track=tC2.begin(); track!=tC2.end(); ++track){
     trackingRecHit_iterator jtB = track->recHitsBegin();
     trackingRecHit_iterator jtE = track->recHitsEnd();
     for (trackingRecHit_iterator jt = jtB;  jt != jtE; ++jt) { 
       const TrackingRecHit* hit = &(**jt);
       rh2[track].push_back(hit);
     }
   }

   if ( (0<tC1.size())&&(0<tC2.size()) ){
    i=-1;
    for (reco::TrackCollection::const_iterator track=tC1.begin(); track!=tC1.end(); ++track){
      i++; 
      if (!selected1[i])continue;
      std::vector<const TrackingRecHit*>& iHits = rh1[track]; 
      unsigned nh1 = iHits.size();
      int qualityMaskT1 = track->qualityMask();
      int j=-1;
      for (reco::TrackCollection::const_iterator track2=tC2.begin(); track2!=tC2.end(); ++track2){
        j++;
        if ((!selected2[j])||(!selected1[i]))continue;
        std::vector<const TrackingRecHit*>& jHits = rh2[track2]; 
        unsigned nh2 = jHits.size();
        int noverlap=0;
        int firstoverlap=0;
        for ( unsigned ih=0; ih<nh1; ++ih ) { 
          const TrackingRecHit* it = iHits[ih];
          if (it->isValid()){
            int jj=-1;
            for ( unsigned jh=0; jh<nh2; ++jh ) { 
              const TrackingRecHit* jt = jHits[jh];
              jj++;
              if (jt->isValid()){
               if (!use_sharesInput){
                float delta = fabs ( it->localPosition().x()-jt->localPosition().x() ); 
                if ((it->geographicalId()==jt->geographicalId())&&(delta<epsilon)) {
                  noverlap++;
                  if ( allowFirstHitShare && ( ih == 0 ) && ( jh == 0 ) ) firstoverlap=1;
                }
               }else{
                if ( it->sharesInput(jt,TrackingRecHit::some) ) {
                  noverlap++;
                  if ( allowFirstHitShare && ( ih == 0 ) && ( jh == 0 ) ) firstoverlap=1;
                }
               }
              }
            }
          }
        }
        int newQualityMask = (qualityMaskT1 | track2->qualityMask()); // take OR of trackQuality 
        int nhit1 = track->numberOfValidHits();
        int nhit2 = track2->numberOfValidHits();
        //std::cout << " trk1 trk2 nhits1 nhits2 nover " << i << " " << j << " " << track->numberOfValidHits() << " "  << track2->numberOfValidHits() << " " << noverlap << " " << fi << " " << fj  <<std::endl;
        if ( (noverlap-firstoverlap) > (std::min(nhit1,nhit2)-firstoverlap)*shareFrac ) {
          if ( nhit1 > nhit2 ){
            selected2[j]=0; 
            selected1[i]=10+newQualityMask; // add 10 to avoid the case where mask = 1
            //std::cout << " removing L2 trk in pair " << std::endl;
          }else{
            if ( nhit1 < nhit2 ){
              selected1[i]=0; 
              selected2[j]=10+newQualityMask;  // add 10 to avoid the case where mask = 1
              //std::cout << " removing L1 trk in pair " << std::endl;
            }else{
              //std::cout << " removing worst chisq in pair " << track->normalizedChi2() << " " << track2->normalizedChi2() << std::endl;
              const double almostSame = 1.001;
              if (track->normalizedChi2() > almostSame * track2->normalizedChi2()) {
                selected1[i]=0;
                selected2[j]=10+newQualityMask; // add 10 to avoid the case where mask = 1
              }else if (track2->normalizedChi2() > almostSame * track->normalizedChi2()) {
                selected2[j]=0;
                selected1[i]=10+newQualityMask; // add 10 to avoid the case where mask = 1
              }else{
                // If tracks from both iterations are virtually identical, choose the one from the first iteration.
                //              std::cout<<"MERGE "<<track->algo()<<" "<<track2->algo()<<" "<<track->normalizedChi2()<<" "<<track2->normalizedChi2()<<" "<<(track->normalizedChi2()-track2->normalizedChi2())/track->normalizedChi2()<<std::endl;
                if (track->algo() <= track2->algo()) {
                  selected2[j]=0;
                  selected1[i]=10+newQualityMask; // add 10 to avoid the case where mask = 1
                }else{
                  selected1[i]=0;
                  selected2[j]=10+newQualityMask; // add 10 to avoid the case where mask = 1
                }
              }
            }//end fi > or = fj
          }//end fi < fj
        }//end got a duplicate
      }//end track2 loop
    }//end track loop
   }//end more than 1 track

  //
  //  output selected tracks - if any
  //
   trackRefs.resize(tC1.size()+tC2.size());
   std::vector<edm::RefToBase<TrajectorySeed> > seedsRefs(tC1.size()+tC2.size());
   size_t current = 0;

   if ( 0<tC1.size() ){
     i=0;
     for (reco::TrackCollection::const_iterator track=tC1.begin(); track!=tC1.end(); 
          ++track, ++current, ++i){
      if (!selected1[i]){
        trackRefs[current] = reco::TrackRef();
        continue;
      }
      const reco::Track & theTrack = * track;
      //fill the TrackCollection
      outputTrks->push_back( reco::Track( theTrack ) );
      if (selected1[i]>1 && promoteQuality){
        outputTrks->back().setQualityMask(selected1[i]-10);
        outputTrks->back().setQuality(qualityToSet);
      }
      if (copyExtras_) {
          //--------NEW----------
          edm::RefToBase<TrajectorySeed> origSeedRef = theTrack.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=clusterProduct(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.
          //--------NEW----------
          // Fill TrackExtra collection
          outputTrkExtras->push_back( reco::TrackExtra( 
                        theTrack.outerPosition(), theTrack.outerMomentum(), theTrack.outerOk(),
                        theTrack.innerPosition(), theTrack.innerMomentum(), theTrack.innerOk(),
                        theTrack.outerStateCovariance(), theTrack.outerDetId(),
                        theTrack.innerStateCovariance(), theTrack.innerDetId(),
                        theTrack.seedDirection(), origSeedRef ) );
          seedsRefs[current]=origSeedRef;
          outputTrks->back().setExtra( reco::TrackExtraRef( refTrkExtras, outputTrkExtras->size() - 1) );
          reco::TrackExtra & tx = outputTrkExtras->back();
          tx.setResiduals(theTrack.residuals());
          // fill TrackingRecHits
          std::vector<const TrackingRecHit*>& iHits = rh1[track]; 
          unsigned nh1 = iHits.size();
          for ( unsigned ih=0; ih<nh1; ++ih ) { 
            const TrackingRecHit* hit = iHits[ih];
            //for( trackingRecHit_iterator hit = itB; hit != itE; ++hit ) {
            outputTrkHits->push_back( hit->clone() );
            tx.add( TrackingRecHitRef( refTrkHits, outputTrkHits->size() - 1) );
          }
      }
      trackRefs[current] = reco::TrackRef(refTrks, outputTrks->size() - 1);
    

    }//end faux loop over tracks
   }//end more than 0 track

   //Fill the trajectories, etc. for 1st collection
   edm::Handle< std::vector<Trajectory> > hTraj1;
   e.getByLabel(trackProducer1, hTraj1);
   edm::Handle< TrajTrackAssociationCollection > hTTAss1;
   e.getByLabel(trackProducer1, hTTAss1);
   refTrajs    = e.getRefBeforePut< std::vector<Trajectory> >();

   if (!hTraj1.failedToGet() && !hTTAss1.failedToGet()){
   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 = 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] );
       }
     }
   }
   }

   short offset = current; //save offset into trackRefs

   if ( 0<tC2.size() ){
    i=0;
    for (reco::TrackCollection::const_iterator track=tC2.begin(); track!=tC2.end();
         ++track, ++current, ++i){
      if (!selected2[i]){
        trackRefs[current] = reco::TrackRef();
        continue;
      }
      const reco::Track & theTrack = * track;
      //fill the TrackCollection
      outputTrks->push_back( reco::Track( theTrack ) );
      if (selected2[i]>1 && promoteQuality){
        outputTrks->back().setQualityMask(selected2[i]-10);
        outputTrks->back().setQuality(qualityToSet);
      }
      if (copyExtras_) {
          //--------NEW----------
          edm::RefToBase<TrajectorySeed> origSeedRef = theTrack.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=clusterProduct(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.
          //--------NEW----------
          // Fill TrackExtra collection
          outputTrkExtras->push_back( reco::TrackExtra( 
                        theTrack.outerPosition(), theTrack.outerMomentum(), theTrack.outerOk(),
                        theTrack.innerPosition(), theTrack.innerMomentum(), theTrack.innerOk(),
                        theTrack.outerStateCovariance(), theTrack.outerDetId(),
                        theTrack.innerStateCovariance(), theTrack.innerDetId(),
                        theTrack.seedDirection(), origSeedRef ) );
          seedsRefs[current]=origSeedRef;
          outputTrks->back().setExtra( reco::TrackExtraRef( refTrkExtras, outputTrkExtras->size() - 1) );
          reco::TrackExtra & tx = outputTrkExtras->back();
          tx.setResiduals(theTrack.residuals());
          // fill TrackingRecHits
          std::vector<const TrackingRecHit*>& jHits = rh2[track]; 
          unsigned nh2 = jHits.size();
          for ( unsigned jh=0; jh<nh2; ++jh ) { 
            const TrackingRecHit* hit = jHits[jh];
            outputTrkHits->push_back( hit->clone() );
            tx.add( TrackingRecHitRef( refTrkHits, outputTrkHits->size() - 1) );
          }
      }
      trackRefs[current] = reco::TrackRef(refTrks, outputTrks->size() - 1);

    }//end faux loop over tracks
   }//end more than 0 track

   //Fill the trajectories, etc. for 2nd collection
   edm::Handle< std::vector<Trajectory> > hTraj2;
   e.getByLabel(trackProducer2, hTraj2);
   edm::Handle< TrajTrackAssociationCollection > hTTAss2;
   e.getByLabel(trackProducer2, hTTAss2);

   if (!hTraj2.failedToGet() && !hTTAss2.failedToGet()){
   for (size_t i = 0, n = hTraj2->size(); i < n; ++i) {
     edm::Ref< std::vector<Trajectory> > trajRef(hTraj2, i);
     TrajTrackAssociationCollection::const_iterator match = hTTAss2->find(trajRef);
     if (match != hTTAss2->end()) {
       const edm::Ref<reco::TrackCollection> &trkRef = match->val; 
       short oldKey = static_cast<short>(trkRef.key()) + offset;
       if (trackRefs[oldKey].isNonnull()) {
           outputTrajs->push_back( Trajectory(*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] ); 
       }
     }
   }}

    e.put(outputTrks);
    if (copyExtras_) {
        e.put(outputTrkExtras);
        e.put(outputTrkHits);
        if (makeReKeyedSeeds_)
          e.put(outputSeeds);
    }
    e.put(outputTrajs);
    e.put(outputTTAss);
    return;

  }//end produce
}