TrackerElectrons/src/SimpleTrackListMergerGsf.cc

// $Id: SimpleTrackListMergerGsf.cc,v 1.2 2009/03/20 17:13:34 loizides Exp $

#include <memory>
#include <string>
#include <iostream>
#include <cmath>
#include <vector>
#include "MitEdm/TrackerElectrons/interface/SimpleTrackListMergerGsf.h"
#include "DataFormats/TrackerRecHit2D/interface/SiStripMatchedRecHit2DCollection.h"
#include "DataFormats/TrackerRecHit2D/interface/SiStripRecHit2DCollection.h"
#include "DataFormats/TrajectorySeed/interface/TrajectorySeedCollection.h"
#include "DataFormats/TrackCandidate/interface/TrackCandidateCollection.h"
#include "FWCore/Framework/interface/ESHandle.h"
#include "FWCore/Framework/interface/MakerMacros.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"

using namespace mitedm;
using namespace edm;

//--------------------------------------------------------------------------------------------------
SimpleTrackListMergerGsf::SimpleTrackListMergerGsf(edm::ParameterSet const &conf) : 
  conf_(conf)
{
  produces<reco::GsfTrackCollection>();
  produces<reco::TrackExtraCollection>();
  produces<reco::GsfTrackExtraCollection>();
  produces<TrackingRecHitCollection>();
  produces<std::vector<Trajectory> >();
}

//--------------------------------------------------------------------------------------------------
SimpleTrackListMergerGsf::~SimpleTrackListMergerGsf() 
{ 
  // Virtual destructor needed.
}  

//--------------------------------------------------------------------------------------------------
void SimpleTrackListMergerGsf::produce(edm::Event& e, const edm::EventSetup& es)
{
  // Functions that gets called by framework every event.

  // retrieve producer name of input GsfTrackCollection(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");
  
  // 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);

  // get Inputs 
  // if 1 input list doesn't exist, make an empty list, issue a warning, and continue
  // this allows SimpleTrackListMergerGsf 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::GsfTrackCollection *TC1 = 0;
  static const reco::GsfTrackCollection s_empty1, s_empty2;
  edm::Handle<reco::GsfTrackCollection> trackCollection1;
  e.getByLabel(trackProducer1, trackCollection1);
  if (trackCollection1.isValid()) {
    TC1 = trackCollection1.product();
    if (0)
      std::cout << "1st collection " << trackProducer1 
                << " has "<< TC1->size() << " tracks" << std::endl;
  } else {
    TC1 = &s_empty1;
    edm::LogWarning("SimpleTrackListMergerGsf") 
      << "1st GsfTrackCollection " << trackProducer1 
      << " not found; will only clean 2nd GsfTrackCollection " << trackProducer2;
  }
  const reco::GsfTrackCollection tC1 = *TC1;

  const reco::GsfTrackCollection *TC2 = 0;
  edm::Handle<reco::GsfTrackCollection> trackCollection2;
  e.getByLabel(trackProducer2, trackCollection2);
  if (trackCollection2.isValid()) {
    TC2 = trackCollection2.product();
    if (0)
      std::cout << "2nd collection " << trackProducer2 << " has "
                << TC2->size() << " tracks" << std::endl;
  } else {
    TC2 = &s_empty2;
    edm::LogWarning("SimpleTrackListMergerGsf") 
      << "2nd GsfTrackCollection " << trackProducer2 
      << " not found; will only clean 1st GsfTrackCollection " << trackProducer1;
  }
  const reco::GsfTrackCollection tC2 = *TC2;
  
  // Step B: create empty output collection
  outputTrks_ = std::auto_ptr<reco::GsfTrackCollection>(new reco::GsfTrackCollection);
  refTrks_    = e.getRefBeforePut<reco::GsfTrackCollection>();      
  
  outputTrkExtras_    = std::auto_ptr<reco::TrackExtraCollection>(new reco::TrackExtraCollection);
  outputGsfTrkExtras_ = 
    std::auto_ptr<reco::GsfTrackExtraCollection>(new reco::GsfTrackExtraCollection);
  refTrkExtras_       = e.getRefBeforePut<reco::TrackExtraCollection>();
  refGsfTrkExtras_    = e.getRefBeforePut<reco::GsfTrackExtraCollection>();
  outputTrkHits_      = std::auto_ptr<TrackingRecHitCollection>(new TrackingRecHitCollection);
  refTrkHits_         = e.getRefBeforePut<TrackingRecHitCollection>();
  outputTrajs_        =  std::auto_ptr< std::vector<Trajectory> >(new std::vector<Trajectory>()); 

  //
  //  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::GsfTrackCollection::const_iterator track=tC1.begin(); track!=tC1.end(); track++) {
      i++;
      if ((short unsigned)track->ndof() < 1){
        selected1[i]=0; 
        if (0)
          std::cout << "L1Track "<< i 
                    << " rejected in SimpleTrackListMergerGsf; ndof() < 1" << std::endl;
        continue;
      }
      if (track->normalizedChi2() > maxNormalizedChisq){
        selected1[i]=0; 
        if (0)
          std::cout << "L1Track "<< i << " rejected in SimpleTrackListMergerGsf; " 
                    << "normalizedChi2() > maxNormalizedChisq " << track->normalizedChi2() 
                    << " " << maxNormalizedChisq << std::endl;
        continue;
      }
      if (track->found() < minFound){
        selected1[i]=0; 
        if (0)
          std::cout << "L1Track "<< i << " rejected in SimpleTrackListMergerGsf; " 
                    << " found() < minFound " << track->found() << " " << minFound << std::endl;
        continue;
      }
      if (track->pt() < minPT){
        selected1[i]=0; 
        if (0)
          std::cout << "L1Track "<< i << " rejected in SimpleTrackListMergerGsf; 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::GsfTrackCollection::const_iterator track=tC2.begin(); track!=tC2.end(); track++) {
      i++;
      if ((short unsigned)track->ndof() < 1){
        selected2[i]=0; 
        if (0) 
        std::cout << "L2Track "<< i << " rejected in SimpleTrackListMergerGsf; ndof() < 1" 
                  << std::endl;
        continue;
      }
      if (track->normalizedChi2() > maxNormalizedChisq){
        selected2[i]=0; 
        if (0)
          std::cout << "L2Track "<< i << " rejected in SimpleTrackListMergerGsf; " 
                    << "normalizedChi2() > maxNormalizedChisq " << track->normalizedChi2() 
                    << " " << maxNormalizedChisq << std::endl;
        continue;
      }
      if (track->found() < minFound){
        selected2[i]=0; 
        if (0)
          std::cout << "L2Track "<< i << " rejected in SimpleTrackListMergerGsf; " 
                    << "found() < minFound " << track->found() << " " << minFound << std::endl;
        continue;
      }
      if (track->pt() < minPT){
        selected2[i]=0; 
        if (0) 
          std::cout << "L2Track "<< i << " rejected in SimpleTrackListMergerGsf; " 
                    << "pt() < minPT " << track->pt() << " " << minPT << std::endl;
        continue;
      }
    } //end loop over tracks
  } //end more than 0 track

  //DON'T BELIEVE WE NEED THIS LOOP
  if (0) {
    //
    //  L1 has > 1 track - try merging
    //
    if (1<tC1.size()) {
      i=-1;
      for (reco::GsfTrackCollection::const_iterator track=tC1.begin(); track!=tC1.end(); track++) {
        i++; 
        if (0)
          std::cout << "Track number "<< i << std::endl; 
        if (!selected1[i]) 
          continue;
        int j=-1;
        for (reco::GsfTrackCollection::const_iterator track2=tC1.begin(); 
             track2!=tC1.end(); track2++) {
          j++;
          if ((j<=i)||(!selected1[j])||(!selected1[i])) 
            continue;
          int noverlap=0;
          for (trackingRecHit_iterator it = track->recHitsBegin();  
               it != track->recHitsEnd(); it++) {
            if ((*it)->isValid()) {
              for (trackingRecHit_iterator jt = track2->recHitsBegin();  
                   jt != track2->recHitsEnd(); jt++){
                if ((*jt)->isValid()) {
                  if (!use_sharesInput){
                    float delta = fabs ((*it)->localPosition().x()-(*jt)->localPosition().x()); 
                    if (((*it)->geographicalId()==(*jt)->geographicalId())&&(delta<epsilon))
                      noverlap++;
                  } else {
                    const TrackingRecHit* kt = &(**jt);
                    if ((*it)->sharesInput(kt,TrackingRecHit::some))
                      noverlap++;
                  }
                }
              }
            }
          }
          float fi=float(noverlap)/float(track->numberOfValidHits()); 
          float fj=float(noverlap)/float(track2->numberOfValidHits());
          if (0)
            std::cout << " trk1 trk2 nhits1 nhits2 nover " << i << " " << j << " " 
                      << track->numberOfValidHits() << " "  << track2->numberOfValidHits() 
                      << " " << noverlap << " " << fi << " " << fj  <<std::endl;
          if ((fi>shareFrac)||(fj>shareFrac)) {
            if (fi<fj) {
              selected1[j]=0;
              if (0)
                std::cout << " removing 2nd trk in pair " << std::endl;
            } else {
              if (fi>fj) {
                selected1[i]=0; 
                if (0) 
                  std::cout << " removing 1st trk in pair " << std::endl;
              } else {
                if (0)
                  std::cout << " removing worst chisq in pair " << std::endl;
                if (track->normalizedChi2() > track2->normalizedChi2()) {
                  selected1[i]=0;
                } else { 
                  selected1[j]=0;
                }
              } //end fi > or = fj
            } //end fi < fj
          } //end got a duplicate
        } //end track2 loop
      } //end track loop
    } //end more than 1 track
  }

  //DON'T BELIEVE WE NEED THIS LOOP EITHER
  if (0) {
    //
    //  L2 has > 1 track - try merging
    //
    if (1<tC2.size()) {
      int i=-1;
      for (reco::GsfTrackCollection::const_iterator track=tC2.begin(); track!=tC2.end(); track++) {
        i++; 
        if (0)
          std::cout << "Track number "<< i << std::endl; 
        if (!selected2[i])
          continue;
        int j=-1;
        for (reco::GsfTrackCollection::const_iterator track2=tC2.begin(); 
             track2!=tC2.end(); track2++){
          j++;
          if ((j<=i)||(!selected2[j])||(!selected2[i])) 
            continue;
          int noverlap=0;
          for (trackingRecHit_iterator it = track->recHitsBegin();
               it != track->recHitsEnd(); it++){
            if ((*it)->isValid()) {
              for (trackingRecHit_iterator jt = track2->recHitsBegin();  
                   jt != track2->recHitsEnd(); jt++) {
                if ((*jt)->isValid()) {
                  if (!use_sharesInput) {
                    float delta = fabs ((*it)->localPosition().x()-(*jt)->localPosition().x()); 
                    if (((*it)->geographicalId()==(*jt)->geographicalId())&&(delta<epsilon))
                      noverlap++;
                  } else {
                    const TrackingRecHit* kt = &(**jt);
                    if ((*it)->sharesInput(kt,TrackingRecHit::some))
                      noverlap++;
                  }
                }
              }
            }
          }
          float fi=float(noverlap)/float(track->numberOfValidHits()); 
          float fj=float(noverlap)/float(track2->numberOfValidHits());
          if (0) 
            std::cout << " trk1 trk2 nhits1 nhits2 nover " << i << " " << j << " " 
                      << track->numberOfValidHits() << " "  << track2->numberOfValidHits() 
                      << " " << noverlap << " " << fi << " " << fj  <<std::endl;
          if ((fi>shareFrac)||(fj>shareFrac)) {
            if (fi<fj) {
              selected2[j]=0; 
              if (0)
                std::cout << " removing 2nd trk in pair " << std::endl;
            } else {
              if (fi>fj) {
                selected2[i]=0; 
                if (0)
                  std::cout << " removing 1st trk in pair " << std::endl;
              } else {
                if (0)
                  std::cout << " removing worst chisq in pair " << std::endl;
                if (track->normalizedChi2() > track2->normalizedChi2()) {
                  selected2[i]=0;
                } else { 
                  selected2[j]=0;
                }
              } //end fi > or = fj
            } //end fi < fj
          } //end got a duplicate
        } //end track2 loop
      } //end track loop
    } //end more than 1 track
  }

  //
  //  L1 and L2 both have > 0 track - try merging
  //
  std::map<reco::GsfTrackCollection::const_iterator, std::vector<const TrackingRecHit*> > rh1;
  std::map<reco::GsfTrackCollection::const_iterator, std::vector<const TrackingRecHit*> > rh2;
  for (reco::GsfTrackCollection::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::GsfTrackCollection::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::GsfTrackCollection::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::GsfTrackCollection::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;
        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++;
                } else {
                  const TrackingRecHit* kt = jt;
                if (it->sharesInput(kt,TrackingRecHit::some))
                  noverlap++;
                }
              }
            }
          }
        }
        int newQualityMask = (qualityMaskT1 | track2->qualityMask()); // take OR of trackQuality 
        float fi=float(noverlap)/float(track->numberOfValidHits()); 
        float fj=float(noverlap)/float(track2->numberOfValidHits());
        if (0)
          std::cout << " trk1 trk2 nhits1 nhits2 nover " << i << " " << j << " " 
                    << track->numberOfValidHits() << " "  << track2->numberOfValidHits() 
                    << " " << noverlap << " " << fi << " " << fj  <<std::endl;
        if ((fi>shareFrac)||(fj>shareFrac)) {
          if (fi<fj){
            selected2[j]=0; 
            selected1[i]=10+newQualityMask; // add 10 to avoid the case where mask = 1
            if (0)
              std::cout << " removing L2 trk in pair " << std::endl;
          } else {
            if (fi>fj) {
              selected1[i]=0; 
              selected2[j]=10+newQualityMask;  // add 10 to avoid the case where mask = 1
              if (0)
                std::cout << " removing L1 trk in pair " << std::endl;
            } else {
              if (0)
                std::cout << " removing worst chisq in pair " << track->normalizedChi2() 
                          << " " << track2->normalizedChi2() << std::endl;
              if (track->normalizedChi2() > track2->normalizedChi2()) {
                selected1[i]=0;
                selected2[j]=10+newQualityMask; // add 10 to avoid the case where mask = 1
              } else {
                selected2[j]=0;
                selected1[i]=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());
  size_t current = 0;
  
  if (0<tC1.size()) {
    i=0;
    for (reco::GsfTrackCollection::const_iterator track=tC1.begin(); track!=tC1.end(); 
         ++track, ++current, ++i){
      if (!selected1[i]) {
        trackRefs_[current] = reco::GsfTrackRef();
        continue;
      }
      const reco::GsfTrack & theTrack = * track;
      //fill the GsfTrackCollection
      outputTrks_->push_back(reco::GsfTrack(theTrack));
      if (selected1[i]>1 && promoteQuality){
        outputTrks_->back().setQualityMask(selected1[i]-10);
        outputTrks_->back().setQuality(qualityToSet);
      }
      // 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(), theTrack.seedRef()));
      outputTrks_->back().setExtra(reco::TrackExtraRef(refTrkExtras_, outputTrkExtras_->size() - 1));

      // Fill GsfTrackExtra collection
      outputGsfTrkExtras_->push_back(*theTrack.gsfExtra());
      outputTrks_->back().setGsfExtra(reco::GsfTrackExtraRef(refGsfTrkExtras_, 
                                                             outputGsfTrkExtras_->size() - 1));

      reco::TrackExtra & tx = outputTrkExtras_->back();
      // 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::GsfTrackRef(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);
 
  outputTrajs_ = std::auto_ptr< std::vector<Trajectory> >(new std::vector<Trajectory>()); 
  refTrajs_    = e.getRefBeforePut< std::vector<Trajectory> >();

  for (size_t i = 0, n = hTraj1->size(); i < n; ++i) {
    edm::Ref< std::vector<Trajectory> > trajRef(hTraj1, i);
    outputTrajs_->push_back(Trajectory(*trajRef));
  }

  if (0<tC2.size()) {
    i=0;
    for (reco::GsfTrackCollection::const_iterator track=tC2.begin(); track!=tC2.end();
         ++track, ++current, ++i){
      if (!selected2[i]){
        trackRefs_[current] = reco::GsfTrackRef();
        continue;
      }
      const reco::GsfTrack & theTrack = * track;
      //fill the GsfTrackCollection
      outputTrks_->push_back(reco::GsfTrack(theTrack));
      if (selected2[i]>1 && promoteQuality){
        outputTrks_->back().setQualityMask(selected2[i]-10);
        outputTrks_->back().setQuality(qualityToSet);
      }
      // 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(), theTrack.seedRef()));
      outputTrks_->back().setExtra(reco::TrackExtraRef(refTrkExtras_, outputTrkExtras_->size() - 1));

      // Fill GsfTrackExtra collection
      outputGsfTrkExtras_->push_back(*theTrack.gsfExtra());
      outputTrks_->back().setGsfExtra(reco::GsfTrackExtraRef(refGsfTrkExtras_, 
                                                            outputGsfTrkExtras_->size() - 1));
      
      reco::TrackExtra &tx = outputTrkExtras_->back();
      // 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::GsfTrackRef(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);

  for (size_t i = 0, n = hTraj2->size(); i < n; ++i) {
    edm::Ref< std::vector<Trajectory> > trajRef(hTraj2, i);
    outputTrajs_->push_back(Trajectory(*trajRef));
  }
   
  e.put(outputTrks_);
  e.put(outputTrkExtras_);
  e.put(outputGsfTrkExtras_);
  e.put(outputTrkHits_);
  e.put(outputTrajs_);
} //end produce

// Define this as a plug-in
DEFINE_FWK_MODULE(SimpleTrackListMergerGsf);
Revision:	1.3
Committed:	Mon Apr 6 19:36:01 2009 UTC (16 years, 1 month ago) by loizides
Content type:	text/plain
Branch:	MAIN
CVS Tags:	Mit_032, Mit_031, Mit_025c_branch2, Mit_025c_branch1, Mit_030, Mit_029c, Mit_029b, Mit_030_pre1, Mit_029a, Mit_029, Mit_029_pre1, Mit_028a, Mit_025c_branch0, Mit_028, Mit_027a, Mit_027, Mit_026, Mit_025e, Mit_025d, Mit_025c, Mit_025b, Mit_025a, Mit_025, Mit_025pre2, Mit_024b, Mit_025pre1, Mit_024a, Mit_024, Mit_023, Mit_022a, Mit_022, Mit_020d, TMit_020d, Mit_020c, Mit_021, Mit_021pre2, Mit_021pre1, Mit_020b, Mit_020a, Mit_020, Mit_020pre1, Mit_018, Mit_017, Mit_017pre3, Mit_017pre2, Mit_017pre1, V07-05-00, Mit_016, Mit_015b, Mit_015a, Mit_015, Mit_014e, Mit_014d, Mit_014c, Mit_014b, ConvRejection-10-06-09, Mit_014a, Mit_014, Mit_014pre3, Mit_014pre2, Mit_014pre1, Mit_013d, Mit_013c, Mit_013b, Mit_013a, Mit_013, Mit_013pre1, Mit_012i, Mit_012h, Mit_012g, Mit_012f, Mit_012e, Mit_012d, Mit_012c, Mit_012b, Mit_012a, Mit_012, Mit_011a, Mit_011, Mit_010a, Mit_010, Mit_009c, Mit_009b, Mit_009a, Mit_009, HEAD
Branch point for:	Mit_025c_branch
Changes since 1.2:	+56 -56 lines
Error occurred while calculating annotation data.
Log Message:	Coding conv.
#	Content
1	// $Id: SimpleTrackListMergerGsf.cc,v 1.2 2009/03/20 17:13:34 loizides Exp $
2
3	#include <memory>
4	#include <string>
5	#include <iostream>
6	#include <cmath>
7	#include <vector>
8	#include "MitEdm/TrackerElectrons/interface/SimpleTrackListMergerGsf.h"
9	#include "DataFormats/TrackerRecHit2D/interface/SiStripMatchedRecHit2DCollection.h"
10	#include "DataFormats/TrackerRecHit2D/interface/SiStripRecHit2DCollection.h"
11	#include "DataFormats/TrajectorySeed/interface/TrajectorySeedCollection.h"
12	#include "DataFormats/TrackCandidate/interface/TrackCandidateCollection.h"
13	#include "FWCore/Framework/interface/ESHandle.h"
14	#include "FWCore/Framework/interface/MakerMacros.h"
15	#include "FWCore/MessageLogger/interface/MessageLogger.h"
16	#include "Geometry/CommonDetUnit/interface/GeomDetUnit.h"
17	#include "Geometry/CommonDetUnit/interface/GeomDet.h"
18	#include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
19	#include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
20
21	using namespace mitedm;
22	using namespace edm;
23
24	//--------------------------------------------------------------------------------------------------
25	SimpleTrackListMergerGsf::SimpleTrackListMergerGsf(edm::ParameterSet const &conf) :
26	conf_(conf)
27	{
28	produces<reco::GsfTrackCollection>();
29	produces<reco::TrackExtraCollection>();
30	produces<reco::GsfTrackExtraCollection>();
31	produces<TrackingRecHitCollection>();
32	produces<std::vector<Trajectory> >();
33	}
34
35	//--------------------------------------------------------------------------------------------------
36	SimpleTrackListMergerGsf::~SimpleTrackListMergerGsf()
37	{
38	// Virtual destructor needed.
39	}
40
41	//--------------------------------------------------------------------------------------------------
42	void SimpleTrackListMergerGsf::produce(edm::Event& e, const edm::EventSetup& es)
43	{
44	// Functions that gets called by framework every event.
45
46	// retrieve producer name of input GsfTrackCollection(s)
47	std::string trackProducer1 = conf_.getParameter<std::string>("TrackProducer1");
48	std::string trackProducer2 = conf_.getParameter<std::string>("TrackProducer2");
49
50	double maxNormalizedChisq = conf_.getParameter<double>("MaxNormalizedChisq");
51	double minPT = conf_.getParameter<double>("MinPT");
52	unsigned int minFound = (unsigned int)conf_.getParameter<int>("MinFound");
53	double epsilon = conf_.getParameter<double>("Epsilon");
54	bool use_sharesInput = true;
55	if (epsilon > 0.0)
56	use_sharesInput = false;
57	double shareFrac = conf_.getParameter<double>("ShareFrac");
58
59	bool promoteQuality = conf_.getParameter<bool>("promoteTrackQuality");
60
61	// New track quality should be read from the file
62	std::string qualityStr = conf_.getParameter<std::string>("newQuality");
63	reco::TrackBase::TrackQuality qualityToSet;
64	if (qualityStr != "") {
65	qualityToSet = reco::TrackBase::qualityByName(conf_.getParameter<std::string>("newQuality"));
66	}
67	else
68	qualityToSet = reco::TrackBase::undefQuality;
69
70	// extract tracker geometry
71	//
72	edm::ESHandle<TrackerGeometry> theG;
73	es.get<TrackerDigiGeometryRecord>().get(theG);
74
75	// get Inputs
76	// if 1 input list doesn't exist, make an empty list, issue a warning, and continue
77	// this allows SimpleTrackListMergerGsf to be used as a cleaner only if handed just one list
78	// if both input lists don't exist, will issue 2 warnings and generate an empty output
79	// collection
80	const reco::GsfTrackCollection *TC1 = 0;
81	static const reco::GsfTrackCollection s_empty1, s_empty2;
82	edm::Handle<reco::GsfTrackCollection> trackCollection1;
83	e.getByLabel(trackProducer1, trackCollection1);
84	if (trackCollection1.isValid()) {
85	TC1 = trackCollection1.product();
86	if (0)
87	std::cout << "1st collection " << trackProducer1
88	<< " has "<< TC1->size() << " tracks" << std::endl;
89	} else {
90	TC1 = &s_empty1;
91	edm::LogWarning("SimpleTrackListMergerGsf")
92	<< "1st GsfTrackCollection " << trackProducer1
93	<< " not found; will only clean 2nd GsfTrackCollection " << trackProducer2;
94	}
95	const reco::GsfTrackCollection tC1 = *TC1;
96
97	const reco::GsfTrackCollection *TC2 = 0;
98	edm::Handle<reco::GsfTrackCollection> trackCollection2;
99	e.getByLabel(trackProducer2, trackCollection2);
100	if (trackCollection2.isValid()) {
101	TC2 = trackCollection2.product();
102	if (0)
103	std::cout << "2nd collection " << trackProducer2 << " has "
104	<< TC2->size() << " tracks" << std::endl;
105	} else {
106	TC2 = &s_empty2;
107	edm::LogWarning("SimpleTrackListMergerGsf")
108	<< "2nd GsfTrackCollection " << trackProducer2
109	<< " not found; will only clean 1st GsfTrackCollection " << trackProducer1;
110	}
111	const reco::GsfTrackCollection tC2 = *TC2;
112
113	// Step B: create empty output collection
114	outputTrks_ = std::auto_ptr<reco::GsfTrackCollection>(new reco::GsfTrackCollection);
115	refTrks_ = e.getRefBeforePut<reco::GsfTrackCollection>();
116
117	outputTrkExtras_ = std::auto_ptr<reco::TrackExtraCollection>(new reco::TrackExtraCollection);
118	outputGsfTrkExtras_ =
119	std::auto_ptr<reco::GsfTrackExtraCollection>(new reco::GsfTrackExtraCollection);
120	refTrkExtras_ = e.getRefBeforePut<reco::TrackExtraCollection>();
121	refGsfTrkExtras_ = e.getRefBeforePut<reco::GsfTrackExtraCollection>();
122	outputTrkHits_ = std::auto_ptr<TrackingRecHitCollection>(new TrackingRecHitCollection);
123	refTrkHits_ = e.getRefBeforePut<TrackingRecHitCollection>();
124	outputTrajs_ = std::auto_ptr< std::vector<Trajectory> >(new std::vector<Trajectory>());
125
126	//
127	// quality cuts first
128	//
129	int i;
130
131	std::vector<int> selected1;
132	for (unsigned int i=0; i<tC1.size(); ++i)
133	selected1.push_back(1);
134
135	if (0<tC1.size()) {
136	i=-1;
137	for (reco::GsfTrackCollection::const_iterator track=tC1.begin(); track!=tC1.end(); track++) {
138	i++;
139	if ((short unsigned)track->ndof() < 1){
140	selected1[i]=0;
141	if (0)
142	std::cout << "L1Track "<< i
143	<< " rejected in SimpleTrackListMergerGsf; ndof() < 1" << std::endl;
144	continue;
145	}
146	if (track->normalizedChi2() > maxNormalizedChisq){
147	selected1[i]=0;
148	if (0)
149	std::cout << "L1Track "<< i << " rejected in SimpleTrackListMergerGsf; "
150	<< "normalizedChi2() > maxNormalizedChisq " << track->normalizedChi2()
151	<< " " << maxNormalizedChisq << std::endl;
152	continue;
153	}
154	if (track->found() < minFound){
155	selected1[i]=0;
156	if (0)
157	std::cout << "L1Track "<< i << " rejected in SimpleTrackListMergerGsf; "
158	<< " found() < minFound " << track->found() << " " << minFound << std::endl;
159	continue;
160	}
161	if (track->pt() < minPT){
162	selected1[i]=0;
163	if (0)
164	std::cout << "L1Track "<< i << " rejected in SimpleTrackListMergerGsf; pt() < minPT "
165	<< track->pt() << " " << minPT << std::endl;
166	continue;
167	}
168	} //end loop over tracks
169	} //end more than 0 track
170
171	std::vector<int> selected2;
172	for (unsigned int i=0; i<tC2.size(); ++i)
173	selected2.push_back(1);
174
175	if (0<tC2.size()) {
176	i=-1;
177	for (reco::GsfTrackCollection::const_iterator track=tC2.begin(); track!=tC2.end(); track++) {
178	i++;
179	if ((short unsigned)track->ndof() < 1){
180	selected2[i]=0;
181	if (0)
182	std::cout << "L2Track "<< i << " rejected in SimpleTrackListMergerGsf; ndof() < 1"
183	<< std::endl;
184	continue;
185	}
186	if (track->normalizedChi2() > maxNormalizedChisq){
187	selected2[i]=0;
188	if (0)
189	std::cout << "L2Track "<< i << " rejected in SimpleTrackListMergerGsf; "
190	<< "normalizedChi2() > maxNormalizedChisq " << track->normalizedChi2()
191	<< " " << maxNormalizedChisq << std::endl;
192	continue;
193	}
194	if (track->found() < minFound){
195	selected2[i]=0;
196	if (0)
197	std::cout << "L2Track "<< i << " rejected in SimpleTrackListMergerGsf; "
198	<< "found() < minFound " << track->found() << " " << minFound << std::endl;
199	continue;
200	}
201	if (track->pt() < minPT){
202	selected2[i]=0;
203	if (0)
204	std::cout << "L2Track "<< i << " rejected in SimpleTrackListMergerGsf; "
205	<< "pt() < minPT " << track->pt() << " " << minPT << std::endl;
206	continue;
207	}
208	} //end loop over tracks
209	} //end more than 0 track
210
211	//DON'T BELIEVE WE NEED THIS LOOP
212	if (0) {
213	//
214	// L1 has > 1 track - try merging
215	//
216	if (1<tC1.size()) {
217	i=-1;
218	for (reco::GsfTrackCollection::const_iterator track=tC1.begin(); track!=tC1.end(); track++) {
219	i++;
220	if (0)
221	std::cout << "Track number "<< i << std::endl;
222	if (!selected1[i])
223	continue;
224	int j=-1;
225	for (reco::GsfTrackCollection::const_iterator track2=tC1.begin();
226	track2!=tC1.end(); track2++) {
227	j++;
228	if ((j<=i)\|\|(!selected1[j])\|\|(!selected1[i]))
229	continue;
230	int noverlap=0;
231	for (trackingRecHit_iterator it = track->recHitsBegin();
232	it != track->recHitsEnd(); it++) {
233	if ((*it)->isValid()) {
234	for (trackingRecHit_iterator jt = track2->recHitsBegin();
235	jt != track2->recHitsEnd(); jt++){
236	if ((*jt)->isValid()) {
237	if (!use_sharesInput){
238	float delta = fabs ((it)->localPosition().x()-(jt)->localPosition().x());
239	if (((it)->geographicalId()==(jt)->geographicalId())&&(delta<epsilon))
240	noverlap++;
241	} else {
242	const TrackingRecHit* kt = &(**jt);
243	if ((*it)->sharesInput(kt,TrackingRecHit::some))
244	noverlap++;
245	}
246	}
247	}
248	}
249	}
250	float fi=float(noverlap)/float(track->numberOfValidHits());
251	float fj=float(noverlap)/float(track2->numberOfValidHits());
252	if (0)
253	std::cout << " trk1 trk2 nhits1 nhits2 nover " << i << " " << j << " "
254	<< track->numberOfValidHits() << " " << track2->numberOfValidHits()
255	<< " " << noverlap << " " << fi << " " << fj <<std::endl;
256	if ((fi>shareFrac)\|\|(fj>shareFrac)) {
257	if (fi<fj) {
258	selected1[j]=0;
259	if (0)
260	std::cout << " removing 2nd trk in pair " << std::endl;
261	} else {
262	if (fi>fj) {
263	selected1[i]=0;
264	if (0)
265	std::cout << " removing 1st trk in pair " << std::endl;
266	} else {
267	if (0)
268	std::cout << " removing worst chisq in pair " << std::endl;
269	if (track->normalizedChi2() > track2->normalizedChi2()) {
270	selected1[i]=0;
271	} else {
272	selected1[j]=0;
273	}
274	} //end fi > or = fj
275	} //end fi < fj
276	} //end got a duplicate
277	} //end track2 loop
278	} //end track loop
279	} //end more than 1 track
280	}
281
282	//DON'T BELIEVE WE NEED THIS LOOP EITHER
283	if (0) {
284	//
285	// L2 has > 1 track - try merging
286	//
287	if (1<tC2.size()) {
288	int i=-1;
289	for (reco::GsfTrackCollection::const_iterator track=tC2.begin(); track!=tC2.end(); track++) {
290	i++;
291	if (0)
292	std::cout << "Track number "<< i << std::endl;
293	if (!selected2[i])
294	continue;
295	int j=-1;
296	for (reco::GsfTrackCollection::const_iterator track2=tC2.begin();
297	track2!=tC2.end(); track2++){
298	j++;
299	if ((j<=i)\|\|(!selected2[j])\|\|(!selected2[i]))
300	continue;
301	int noverlap=0;
302	for (trackingRecHit_iterator it = track->recHitsBegin();
303	it != track->recHitsEnd(); it++){
304	if ((*it)->isValid()) {
305	for (trackingRecHit_iterator jt = track2->recHitsBegin();
306	jt != track2->recHitsEnd(); jt++) {
307	if ((*jt)->isValid()) {
308	if (!use_sharesInput) {
309	float delta = fabs ((it)->localPosition().x()-(jt)->localPosition().x());
310	if (((it)->geographicalId()==(jt)->geographicalId())&&(delta<epsilon))
311	noverlap++;
312	} else {
313	const TrackingRecHit* kt = &(**jt);
314	if ((*it)->sharesInput(kt,TrackingRecHit::some))
315	noverlap++;
316	}
317	}
318	}
319	}
320	}
321	float fi=float(noverlap)/float(track->numberOfValidHits());
322	float fj=float(noverlap)/float(track2->numberOfValidHits());
323	if (0)
324	std::cout << " trk1 trk2 nhits1 nhits2 nover " << i << " " << j << " "
325	<< track->numberOfValidHits() << " " << track2->numberOfValidHits()
326	<< " " << noverlap << " " << fi << " " << fj <<std::endl;
327	if ((fi>shareFrac)\|\|(fj>shareFrac)) {
328	if (fi<fj) {
329	selected2[j]=0;
330	if (0)
331	std::cout << " removing 2nd trk in pair " << std::endl;
332	} else {
333	if (fi>fj) {
334	selected2[i]=0;
335	if (0)
336	std::cout << " removing 1st trk in pair " << std::endl;
337	} else {
338	if (0)
339	std::cout << " removing worst chisq in pair " << std::endl;
340	if (track->normalizedChi2() > track2->normalizedChi2()) {
341	selected2[i]=0;
342	} else {
343	selected2[j]=0;
344	}
345	} //end fi > or = fj
346	} //end fi < fj
347	} //end got a duplicate
348	} //end track2 loop
349	} //end track loop
350	} //end more than 1 track
351	}
352
353	//
354	// L1 and L2 both have > 0 track - try merging
355	//
356	std::map<reco::GsfTrackCollection::const_iterator, std::vector<const TrackingRecHit*> > rh1;
357	std::map<reco::GsfTrackCollection::const_iterator, std::vector<const TrackingRecHit*> > rh2;
358	for (reco::GsfTrackCollection::const_iterator track=tC1.begin(); track!=tC1.end(); ++track) {
359	trackingRecHit_iterator itB = track->recHitsBegin();
360	trackingRecHit_iterator itE = track->recHitsEnd();
361	for (trackingRecHit_iterator it = itB; it != itE; ++it) {
362	const TrackingRecHit* hit = &(**it);
363	rh1[track].push_back(hit);
364	}
365	}
366	for (reco::GsfTrackCollection::const_iterator track=tC2.begin(); track!=tC2.end(); ++track){
367	trackingRecHit_iterator jtB = track->recHitsBegin();
368	trackingRecHit_iterator jtE = track->recHitsEnd();
369	for (trackingRecHit_iterator jt = jtB; jt != jtE; ++jt) {
370	const TrackingRecHit* hit = &(**jt);
371	rh2[track].push_back(hit);
372	}
373	}
374
375	if ((0<tC1.size())&&(0<tC2.size())) {
376	i=-1;
377	for (reco::GsfTrackCollection::const_iterator track=tC1.begin(); track!=tC1.end(); ++track) {
378	i++;
379	if (!selected1[i])
380	continue;
381	std::vector<const TrackingRecHit*>& iHits = rh1[track];
382	unsigned nh1 = iHits.size();
383	int qualityMaskT1 = track->qualityMask();
384	int j=-1;
385	for (reco::GsfTrackCollection::const_iterator track2=tC2.begin();
386	track2!=tC2.end(); ++track2) {
387	j++;
388	if ((!selected2[j])\|\|(!selected1[i]))
389	continue;
390	std::vector<const TrackingRecHit*>& jHits = rh2[track2];
391	unsigned nh2 = jHits.size();
392	int noverlap=0;
393	for (unsigned ih=0; ih<nh1; ++ih) {
394	const TrackingRecHit* it = iHits[ih];
395	if (it->isValid()) {
396	int jj=-1;
397	for (unsigned jh=0; jh<nh2; ++jh) {
398	const TrackingRecHit* jt = jHits[jh];
399	jj++;
400	if (jt->isValid()) {
401	if (!use_sharesInput) {
402	float delta = fabs(it->localPosition().x()-jt->localPosition().x());
403	if ((it->geographicalId()==jt->geographicalId())&&(delta<epsilon))
404	noverlap++;
405	} else {
406	const TrackingRecHit* kt = jt;
407	if (it->sharesInput(kt,TrackingRecHit::some))
408	noverlap++;
409	}
410	}
411	}
412	}
413	}
414	int newQualityMask = (qualityMaskT1 \| track2->qualityMask()); // take OR of trackQuality
415	float fi=float(noverlap)/float(track->numberOfValidHits());
416	float fj=float(noverlap)/float(track2->numberOfValidHits());
417	if (0)
418	std::cout << " trk1 trk2 nhits1 nhits2 nover " << i << " " << j << " "
419	<< track->numberOfValidHits() << " " << track2->numberOfValidHits()
420	<< " " << noverlap << " " << fi << " " << fj <<std::endl;
421	if ((fi>shareFrac)\|\|(fj>shareFrac)) {
422	if (fi<fj){
423	selected2[j]=0;
424	selected1[i]=10+newQualityMask; // add 10 to avoid the case where mask = 1
425	if (0)
426	std::cout << " removing L2 trk in pair " << std::endl;
427	} else {
428	if (fi>fj) {
429	selected1[i]=0;
430	selected2[j]=10+newQualityMask; // add 10 to avoid the case where mask = 1
431	if (0)
432	std::cout << " removing L1 trk in pair " << std::endl;
433	} else {
434	if (0)
435	std::cout << " removing worst chisq in pair " << track->normalizedChi2()
436	<< " " << track2->normalizedChi2() << std::endl;
437	if (track->normalizedChi2() > track2->normalizedChi2()) {
438	selected1[i]=0;
439	selected2[j]=10+newQualityMask; // add 10 to avoid the case where mask = 1
440	} else {
441	selected2[j]=0;
442	selected1[i]=10+newQualityMask; // add 10 to avoid the case where mask = 1
443	}
444	} //end fi > or = fj
445	} //end fi < fj
446	} //end got a duplicate
447	} //end track2 loop
448	} //end track loop
449	} //end more than 1 track
450
451	//
452	// output selected tracks - if any
453	//
454	trackRefs_.resize(tC1.size()+tC2.size());
455	size_t current = 0;
456
457	if (0<tC1.size()) {
458	i=0;
459	for (reco::GsfTrackCollection::const_iterator track=tC1.begin(); track!=tC1.end();
460	++track, ++current, ++i){
461	if (!selected1[i]) {
462	trackRefs_[current] = reco::GsfTrackRef();
463	continue;
464	}
465	const reco::GsfTrack & theTrack = * track;
466	//fill the GsfTrackCollection
467	outputTrks_->push_back(reco::GsfTrack(theTrack));
468	if (selected1[i]>1 && promoteQuality){
469	outputTrks_->back().setQualityMask(selected1[i]-10);
470	outputTrks_->back().setQuality(qualityToSet);
471	}
472	// Fill TrackExtra collection
473	outputTrkExtras_->push_back(reco::TrackExtra(
474	theTrack.outerPosition(), theTrack.outerMomentum(), theTrack.outerOk(),
475	theTrack.innerPosition(), theTrack.innerMomentum(), theTrack.innerOk(),
476	theTrack.outerStateCovariance(), theTrack.outerDetId(),
477	theTrack.innerStateCovariance(), theTrack.innerDetId(),
478	theTrack.seedDirection(), theTrack.seedRef()));
479	outputTrks_->back().setExtra(reco::TrackExtraRef(refTrkExtras_, outputTrkExtras_->size() - 1));
480
481	// Fill GsfTrackExtra collection
482	outputGsfTrkExtras_->push_back(*theTrack.gsfExtra());
483	outputTrks_->back().setGsfExtra(reco::GsfTrackExtraRef(refGsfTrkExtras_,
484	outputGsfTrkExtras_->size() - 1));
485
486	reco::TrackExtra & tx = outputTrkExtras_->back();
487	// fill TrackingRecHits
488	std::vector<const TrackingRecHit*>& iHits = rh1[track];
489	unsigned nh1 = iHits.size();
490	for (unsigned ih=0; ih<nh1; ++ih) {
491	const TrackingRecHit* hit = iHits[ih];
492	//for(trackingRecHit_iterator hit = itB; hit != itE; ++hit) {
493	outputTrkHits_->push_back(hit->clone());
494	tx.add(TrackingRecHitRef(refTrkHits_, outputTrkHits_->size() - 1));
495	}
496	trackRefs_[current] = reco::GsfTrackRef(refTrks_, outputTrks_->size() - 1);
497
498	} //end faux loop over tracks
499	} //end more than 0 track
500
501	//Fill the trajectories, etc. for 1st collection
502	edm::Handle< std::vector<Trajectory> > hTraj1;
503	e.getByLabel(trackProducer1, hTraj1);
504
505	outputTrajs_ = std::auto_ptr< std::vector<Trajectory> >(new std::vector<Trajectory>());
506	refTrajs_ = e.getRefBeforePut< std::vector<Trajectory> >();
507
508	for (size_t i = 0, n = hTraj1->size(); i < n; ++i) {
509	edm::Ref< std::vector<Trajectory> > trajRef(hTraj1, i);
510	outputTrajs_->push_back(Trajectory(*trajRef));
511	}
512
513	if (0<tC2.size()) {
514	i=0;
515	for (reco::GsfTrackCollection::const_iterator track=tC2.begin(); track!=tC2.end();
516	++track, ++current, ++i){
517	if (!selected2[i]){
518	trackRefs_[current] = reco::GsfTrackRef();
519	continue;
520	}
521	const reco::GsfTrack & theTrack = * track;
522	//fill the GsfTrackCollection
523	outputTrks_->push_back(reco::GsfTrack(theTrack));
524	if (selected2[i]>1 && promoteQuality){
525	outputTrks_->back().setQualityMask(selected2[i]-10);
526	outputTrks_->back().setQuality(qualityToSet);
527	}
528	// Fill TrackExtra collection
529	outputTrkExtras_->push_back(reco::TrackExtra(
530	theTrack.outerPosition(), theTrack.outerMomentum(), theTrack.outerOk(),
531	theTrack.innerPosition(), theTrack.innerMomentum(), theTrack.innerOk(),
532	theTrack.outerStateCovariance(), theTrack.outerDetId(),
533	theTrack.innerStateCovariance(), theTrack.innerDetId(),
534	theTrack.seedDirection(), theTrack.seedRef()));
535	outputTrks_->back().setExtra(reco::TrackExtraRef(refTrkExtras_, outputTrkExtras_->size() - 1));
536
537	// Fill GsfTrackExtra collection
538	outputGsfTrkExtras_->push_back(*theTrack.gsfExtra());
539	outputTrks_->back().setGsfExtra(reco::GsfTrackExtraRef(refGsfTrkExtras_,
540	outputGsfTrkExtras_->size() - 1));
541
542	reco::TrackExtra &tx = outputTrkExtras_->back();
543	// fill TrackingRecHits
544	std::vector<const TrackingRecHit*> &jHits = rh2[track];
545	unsigned nh2 = jHits.size();
546	for (unsigned jh=0; jh<nh2; ++jh) {
547	const TrackingRecHit* hit = jHits[jh];
548	outputTrkHits_->push_back(hit->clone());
549	tx.add(TrackingRecHitRef(refTrkHits_, outputTrkHits_->size() - 1));
550	}
551	trackRefs_[current] = reco::GsfTrackRef(refTrks_, outputTrks_->size() - 1);
552
553	} //end faux loop over tracks
554	} //end more than 0 track
555
556	// Fill the trajectories, etc. for 2nd collection
557	edm::Handle< std::vector<Trajectory> > hTraj2;
558	e.getByLabel(trackProducer2, hTraj2);
559
560	for (size_t i = 0, n = hTraj2->size(); i < n; ++i) {
561	edm::Ref< std::vector<Trajectory> > trajRef(hTraj2, i);
562	outputTrajs_->push_back(Trajectory(*trajRef));
563	}
564
565	e.put(outputTrks_);
566	e.put(outputTrkExtras_);
567	e.put(outputGsfTrkExtras_);
568	e.put(outputTrkHits_);
569	e.put(outputTrajs_);
570	} //end produce
571
572	// Define this as a plug-in
573	DEFINE_FWK_MODULE(SimpleTrackListMergerGsf);