GoodSeedProducer.cc

Go to the documentation of this file.

// -*- C++ -*-
//
// Package:    PFTracking
// Class:      GoodSeedProducer
// 
// Original Author:  Michele Pioppi
// March 2010. F. Beaudette. Produce PreId information

#include "RecoParticleFlow/PFTracking/plugins/GoodSeedProducer.h"
#include "RecoParticleFlow/PFTracking/interface/PFTrackTransformer.h"
#include "RecoParticleFlow/PFClusterTools/interface/PFResolutionMap.h"

#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "FWCore/ParameterSet/interface/FileInPath.h"
#include "FWCore/Framework/interface/MakerMacros.h"
#include "DataFormats/EgammaReco/interface/ElectronSeed.h"
#include "DataFormats/EgammaReco/interface/ElectronSeedFwd.h"
#include "DataFormats/TrajectorySeed/interface/TrajectorySeedCollection.h"
#include "DataFormats/TrajectorySeed/interface/PropagationDirection.h"
#include "DataFormats/ParticleFlowReco/interface/PreId.h"
#include "TrackingTools/TrackFitters/interface/TrajectoryFitter.h"
#include "TrackingTools/PatternTools/interface/TrajectorySmoother.h"
#include "TrackingTools/Records/interface/TransientRecHitRecord.h"  
#include "TrackingTools/Records/interface/TrackingComponentsRecord.h" 
#include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
#include "FastSimulation/BaseParticlePropagator/interface/BaseParticlePropagator.h"
#include "MagneticField/Engine/interface/MagneticField.h"
#include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"

#include "DataFormats/Math/interface/deltaR.h"

#include <fstream>
#include <string>
#include "TMath.h"
#include "Math/VectorUtil.h"
#include "TMVA/MethodBDT.h"

using namespace edm;
using namespace std;
using namespace reco;

GoodSeedProducer::GoodSeedProducer(const ParameterSet& iConfig, const goodseedhelpers::HeavyObjectCache*):
  pfTransformer_(nullptr),
  conf_(iConfig),
  resMapEtaECAL_(nullptr),
  resMapPhiECAL_(nullptr)
{
  LogInfo("GoodSeedProducer")<<"Electron PreIdentification started  ";
  
  //now do what ever initialization is needed
  std::vector<edm::InputTag> tags =   iConfig.getParameter< vector < InputTag > >("TkColList");
  for(unsigned int i=0;i<tags.size();++i) {
    trajContainers_.push_back(consumes<vector<Trajectory> >(tags[i]));
    tracksContainers_.push_back(consumes<reco::TrackCollection>(tags[i]));
  }
  
  minPt_=iConfig.getParameter<double>("MinPt");
  maxPt_=iConfig.getParameter<double>("MaxPt");
  maxEta_=iConfig.getParameter<double>("MaxEta");

  HcalIsolWindow_                       =iConfig.getParameter<double>("HcalWindow");
  EcalStripSumE_minClusEnergy_ = iConfig.getParameter<double>("EcalStripSumE_minClusEnergy");
  EcalStripSumE_deltaEta_ = iConfig.getParameter<double>("EcalStripSumE_deltaEta");
  EcalStripSumE_deltaPhiOverQ_minValue_ = iConfig.getParameter<double>("EcalStripSumE_deltaPhiOverQ_minValue");
  EcalStripSumE_deltaPhiOverQ_maxValue_ = iConfig.getParameter<double>("EcalStripSumE_deltaPhiOverQ_maxValue");
  minEoverP_= iConfig.getParameter<double>("EOverPLead_minValue");
  maxHoverP_= iConfig.getParameter<double>("HOverPLead_maxValue");
 
  pfCLusTagECLabel_=consumes<reco::PFClusterCollection>(iConfig.getParameter<InputTag>("PFEcalClusterLabel"));

  pfCLusTagHCLabel_=consumes<reco::PFClusterCollection>(iConfig.getParameter<InputTag>("PFHcalClusterLabel"));  

  pfCLusTagPSLabel_=consumes<reco::PFClusterCollection>(iConfig.getParameter<InputTag>("PFPSClusterLabel"));
  
  preidgsf_ = iConfig.getParameter<string>("PreGsfLabel");
  preidckf_ = iConfig.getParameter<string>("PreCkfLabel");
  preidname_= iConfig.getParameter<string>("PreIdLabel");
  
  
  fitterName_ = iConfig.getParameter<string>("Fitter");
  smootherName_ = iConfig.getParameter<string>("Smoother");
  
  
  nHitsInSeed_=iConfig.getParameter<int>("NHitsInSeed");

  clusThreshold_=iConfig.getParameter<double>("ClusterThreshold");
  
  minEp_=iConfig.getParameter<double>("MinEOverP");
  maxEp_=iConfig.getParameter<double>("MaxEOverP");

  //collection to produce
  produceCkfseed_ = iConfig.getUntrackedParameter<bool>("ProduceCkfSeed",false);

  // to disable the electron part (for HI collisions for examples) 
  disablePreId_ = iConfig.getUntrackedParameter<bool>("DisablePreId",false);  

  producePreId_ = iConfig.getUntrackedParameter<bool>("ProducePreId",true);  
  //  if no electron, cannot produce the preid
  if(disablePreId_) 
    producePreId_=false;
  PtThresholdSavePredId_ = iConfig.getUntrackedParameter<double>("PtThresholdSavePreId",1.);  

  LogDebug("GoodSeedProducer")<<"Seeds for GSF will be produced ";

  // no disablePreId_ switch here. The collection will be empty if it is true
  produces<ElectronSeedCollection>(preidgsf_);

  if(produceCkfseed_){
    LogDebug("GoodSeedProducer")<<"Seeds for CKF will be produced ";
    produces<TrajectorySeedCollection>(preidckf_);
  }

  if(producePreId_){
    LogDebug("GoodSeedProducer")<<"PreId debugging information will be produced ";

    produces<PreIdCollection>(preidname_);
    if(tracksContainers_.size()==1) // do not make a value map if more than one input track collection
      produces<edm::ValueMap<reco::PreIdRef> >(preidname_);
  } 
  
  useQuality_   = iConfig.getParameter<bool>("UseQuality");
  trackQuality_=TrackBase::qualityByName(iConfig.getParameter<std::string>("TrackQuality"));

  useTmva_= iConfig.getUntrackedParameter<bool>("UseTMVA",false);

  Min_dr_ = iConfig.getParameter<double>("Min_dr");

  trackerRecHitBuilderName_ = iConfig.getParameter<std::string>("TTRHBuilder");

}


//
// member functions
//

// ------------ method called to produce the data  ------------
void
GoodSeedProducer::produce(Event& iEvent, const EventSetup& iSetup)
{

  LogDebug("GoodSeedProducer")<<"START event: "<<iEvent.id().event()
                  <<" in run "<<iEvent.id().run();
  //Create empty output collections
  auto output_preid = std::make_unique<ElectronSeedCollection>();
  auto output_nopre = std::make_unique<TrajectorySeedCollection>();
  auto output_preidinfo = std::make_unique<PreIdCollection>();
  auto preIdMap_p = std::make_unique<edm::ValueMap<reco::PreIdRef>>();
  edm::ValueMap<reco::PreIdRef>::Filler mapFiller(*preIdMap_p);

  //Tracking Tools
  if(!disablePreId_)
    {
      edm::ESHandle<TrajectoryFitter> aFitter;
      edm::ESHandle<TrajectorySmoother> aSmoother;
      iSetup.get<TrajectoryFitter::Record>().get(fitterName_, aFitter);
      iSetup.get<TrajectoryFitter::Record>().get(smootherName_, aSmoother);
      smoother_.reset(aSmoother->clone());
      fitter_ = aFitter->clone();
      edm::ESHandle<TransientTrackingRecHitBuilder> theTrackerRecHitBuilder;
      iSetup.get<TransientRecHitRecord>().get(trackerRecHitBuilderName_,theTrackerRecHitBuilder);
      hitCloner = static_cast<TkTransientTrackingRecHitBuilder const *>(theTrackerRecHitBuilder.product())->cloner();
      fitter_->setHitCloner(&hitCloner);
      smoother_->setHitCloner(&hitCloner);
    }

  // clear temporary maps
  refMap_.clear();

  //Magnetic Field
  ESHandle<MagneticField> magneticField;
  iSetup.get<IdealMagneticFieldRecord>().get(magneticField);

  //Handle input collections
  //ECAL clusters         
  Handle<PFClusterCollection> theECPfClustCollection;
  iEvent.getByToken(pfCLusTagECLabel_,theECPfClustCollection);
  

 vector<PFCluster const *> basClus;
  for ( auto const & klus : *theECPfClustCollection.product() ) {
    if(klus.correctedEnergy()>clusThreshold_) basClus.push_back(&klus);
  }

  //HCAL clusters
  Handle<PFClusterCollection> theHCPfClustCollection;
  iEvent.getByToken(pfCLusTagHCLabel_,theHCPfClustCollection);
  
  //PS clusters
  Handle<PFClusterCollection> thePSPfClustCollection;
  iEvent.getByToken(pfCLusTagPSLabel_,thePSPfClustCollection);

  //Vector of track collections
  for (unsigned int istr=0; istr<tracksContainers_.size();++istr){
    
    //Track collection
    Handle<TrackCollection> tkRefCollection;
    iEvent.getByToken(tracksContainers_[istr], tkRefCollection);
    const TrackCollection&  Tk=*(tkRefCollection.product());
    
    LogDebug("GoodSeedProducer")<<"Number of tracks in collection "
                                <<"tracksContainers_[" << istr << "] to be analyzed "
                                <<Tk.size();

    //loop over the track collection
    for(unsigned int i=0;i<Tk.size();++i){      
      if (useQuality_ &&
      (!(Tk[i].quality(trackQuality_)))) continue;
      
      reco::PreId myPreId;
      bool GoodPreId=false;

      TrackRef trackRef(tkRefCollection, i);
      math::XYZVectorF tkmom(Tk[i].momentum());
      auto tketa= tkmom.eta();
      auto tkpt = std::sqrt(tkmom.perp2());
      auto const & Seed=(*trackRef->seedRef());
      if(!disablePreId_)
    {
      int ipteta=getBin(Tk[i].eta(),Tk[i].pt());
      int ibin=ipteta*9;
      
      float oPTOB=1.f/std::sqrt(Tk[i].innerMomentum().mag2()); // FIXME the original code was buggy should be outerMomentum...
      //  float chikfred=Tk[i].normalizedChi2();
      float nchi=Tk[i].normalizedChi2();

      int nhitpi=Tk[i].found();
      float EP=0;
      
      // set track info
      myPreId.setTrack(trackRef);
      //CLUSTERS - TRACK matching
      
      auto pfmass=  0.0005;
      auto pfoutenergy=sqrt((pfmass*pfmass)+Tk[i].outerMomentum().Mag2());

      XYZTLorentzVector mom =XYZTLorentzVector(Tk[i].outerMomentum().x(),
                           Tk[i].outerMomentum().y(),
                           Tk[i].outerMomentum().z(),
                           pfoutenergy);
      XYZTLorentzVector pos =   XYZTLorentzVector(Tk[i].outerPosition().x(),
                              Tk[i].outerPosition().y(),
                              Tk[i].outerPosition().z(),
                              0.);

      BaseParticlePropagator theOutParticle( RawParticle(mom,pos),
                    0,0,B_.z());
      theOutParticle.setCharge(Tk[i].charge());
      
      theOutParticle.propagateToEcalEntrance(false);
      

      
      float toteta=1000.f;
      float totphi=1000.f;
      float dr=1000.f;
      float EE=0.f;
      float feta=0.f;
      GlobalPoint ElecTrkEcalPos(0,0,0);

      PFClusterRef clusterRef;
      math::XYZPoint meanShowerSaved;
      if(theOutParticle.getSuccess()!=0){
         ElecTrkEcalPos=GlobalPoint(theOutParticle.vertex().x(),
                            theOutParticle.vertex().y(),
                    theOutParticle.vertex().z()
                                       );

            constexpr float psLim = 2.50746495928f; // std::sinh(1.65f);
            bool isBelowPS= (ElecTrkEcalPos.z()*ElecTrkEcalPos.z()) > (psLim*psLim)*ElecTrkEcalPos.perp2();
        // bool isBelowPS=(std::abs(ElecTrkEcalPos.eta())>1.65f);   
    
        unsigned clusCounter=0;
        float max_ee = 0;
        for(auto aClus : basClus) {

          float tmp_ep=float(aClus->correctedEnergy())*oPTOB;
              if ((tmp_ep<minEp_)|(tmp_ep>maxEp_)) { ++clusCounter; continue;}
        
          double ecalShowerDepth
        = PFCluster::getDepthCorrection(aClus->correctedEnergy(),
                        isBelowPS,
                        false);
          auto mom = theOutParticle.momentum().Vect();
          auto meanShower = ElecTrkEcalPos +
        GlobalVector(mom.x(),mom.y(),mom.z()).unit()*ecalShowerDepth;   
      
          float etarec=meanShower.eta();
          float phirec=meanShower.phi();
         

          float tmp_phi=std::abs(aClus->positionREP().phi()-phirec);
          if (tmp_phi>float(TMath::Pi())) tmp_phi-= float(TMath::TwoPi());
          
          float tmp_dr=std::sqrt(std::pow(tmp_phi,2.f)+
                std::pow(aClus->positionREP().eta()-etarec,2.f));
      
          if (tmp_dr<dr){
        dr=tmp_dr;
        if(dr < Min_dr_){ // find the most closest and energetic ECAL cluster
          if(aClus->correctedEnergy() > max_ee){

            toteta=aClus->positionREP().eta()-etarec;
            totphi=tmp_phi;
            EP=tmp_ep;
            EE=aClus->correctedEnergy();
            feta= aClus->positionREP().eta();
            clusterRef = PFClusterRef(theECPfClustCollection,clusCounter);
            meanShowerSaved = meanShower;
            
          }
        }
          }
              ++clusCounter;
        }
      }
      float trk_ecalDeta_ = fabs(toteta);
      float trk_ecalDphi_ = fabs(totphi);

      //Resolution maps
      auto ecaletares 
        = resMapEtaECAL_->GetBinContent(resMapEtaECAL_->FindBin(feta,EE));
      auto ecalphires 
        = resMapPhiECAL_->GetBinContent(resMapPhiECAL_->FindBin(feta,EE)); 
      
      //geomatrical compatibility
      float chieta=(toteta!=1000.f)? toteta/ecaletares : toteta;
      float chiphi=(totphi!=1000.f)? totphi/ecalphires : totphi;
      float chichi= sqrt(chieta*chieta + chiphi*chiphi);
      
      //Matching criteria
      float eta_cut = thr[ibin+0];
      float phi_cut = thr[ibin+1];
      float ep_cutmin=thr[ibin+2];
      bool GoodMatching= ((trk_ecalDeta_<eta_cut) && (trk_ecalDphi_<phi_cut) && (EP>ep_cutmin) && (nhitpi>10));

      bool EcalMatching=GoodMatching;
      
      if (tkpt>maxPt_) GoodMatching=true;
      if (tkpt<minPt_) GoodMatching=false;


  
      math::XYZPoint myPoint(ElecTrkEcalPos.x(),ElecTrkEcalPos.y(),ElecTrkEcalPos.z());
      myPreId.setECALMatchingProperties(clusterRef,myPoint,meanShowerSaved,std::abs(toteta),std::abs(totphi),chieta,
                        chiphi,chichi,EP);
      myPreId.setECALMatching(EcalMatching);


      bool GoodRange= ((std::abs(tketa)<maxEta_) & 
               (tkpt>minPt_));
      //KF FILTERING FOR UNMATCHED EVENTS
      int hit1max=int(thr[ibin+3]);
      float chiredmin=thr[ibin+4];
      bool GoodKFFiltering =
        ((nchi>chiredmin) | (nhitpi<hit1max));
      

      myPreId.setTrackFiltering(GoodKFFiltering);

      bool GoodTkId= false;
      
      if((!GoodMatching) &&(GoodKFFiltering) &&(GoodRange)){
        chired=1000;
        chiRatio=1000;
        dpt=0;
        nhit=nhitpi;
        chikfred = nchi;
        trk_ecalDeta = trk_ecalDeta_;
        trk_ecalDphi = trk_ecalDphi_;
      
        Trajectory::ConstRecHitContainer tmp;
            auto hb = Tk[i].recHitsBegin();
            for(unsigned int h=0;h<Tk[i].recHitsSize();h++){
              auto recHit = *(hb+h); tmp.push_back(recHit->cloneSH());
            }
            auto const & theTrack = Tk[i]; 
            GlobalVector gv(theTrack.innerMomentum().x(),theTrack.innerMomentum().y(),theTrack.innerMomentum().z());
            GlobalPoint  gp(theTrack.innerPosition().x(),theTrack.innerPosition().y(),theTrack.innerPosition().z());
            GlobalTrajectoryParameters gtps(gp,gv,theTrack.charge(),&*magneticField);
            TrajectoryStateOnSurface tsos(gtps,theTrack.innerStateCovariance(),*tmp[0]->surface());
        Trajectory  && FitTjs= fitter_->fitOne(Seed,tmp,tsos);
    
          if(FitTjs.isValid()){
        Trajectory && SmooTjs= smoother_->trajectory(FitTjs);
          if(SmooTjs.isValid()){
        
            //Track refitted with electron hypothesis
        
            float pt_out=SmooTjs.firstMeasurement().
              updatedState().globalMomentum().perp();
            float pt_in=SmooTjs.lastMeasurement().
              updatedState().globalMomentum().perp();
            dpt=(pt_in>0) ? fabs(pt_out-pt_in)/pt_in : 0.;
            // the following is simply the number of degrees of freedom
            chiRatio=SmooTjs.chiSquared()/Tk[i].chi2();
            chired=chiRatio*chikfred;

          }
        }
         
    
        //TMVA Analysis
        if(useTmva_){
    
          eta=tketa;
          pt=tkpt;
          eP=EP;
              float vars[10] = { nhit, chikfred, dpt, eP, chiRatio, chired, trk_ecalDeta, trk_ecalDphi, pt, eta};
              
          float Ytmva = globalCache()->gbr[ipteta]->GetClassifier( vars );
          
          float BDTcut=thr[ibin+5]; 
          if ( Ytmva>BDTcut) GoodTkId=true;
          myPreId.setMVA(GoodTkId,Ytmva);
          myPreId.setTrackProperties(chired,chiRatio,dpt);
        }else{ 
              
          float chiratiocut=thr[ibin+6]; 
          float gschicut=thr[ibin+7]; 
          float gsptmin=thr[ibin+8];

          GoodTkId=((dpt>gsptmin)&(chired<gschicut)&(chiRatio<chiratiocut));      
       
        }
      }
    
      GoodPreId= GoodTkId | GoodMatching; 

      myPreId.setFinalDecision(GoodPreId);

#ifdef EDM_ML_DEBUG      
      if(GoodPreId)
        LogDebug("GoodSeedProducer")<<"Track (pt= "<<Tk[i].pt()<<
          "GeV/c, eta= "<<Tk[i].eta() <<
          ") preidentified for agreement between  track and ECAL cluster";
      if(GoodPreId &&(!GoodMatching))
        LogDebug("GoodSeedProducer")<<"Track (pt= "<<Tk[i].pt()<<
          "GeV/c, eta= "<<Tk[i].eta() <<
          ") preidentified only for track properties";
#endif

    } // end of !disablePreId_

      
      if (GoodPreId){
    //NEW SEED with n hits  
    ElectronSeed NewSeed(Seed);
    NewSeed.setCtfTrack(trackRef);
    output_preid->push_back(NewSeed);
      }else{
    if (produceCkfseed_){
      output_nopre->push_back(Seed);
    }
      }
      if(producePreId_ && myPreId.pt()>PtThresholdSavePredId_)
    {
      // save the index of the PreId object as to be able to create a Ref later
      refMap_[trackRef] = output_preidinfo->size();
      output_preidinfo->push_back(myPreId);   
    }
    } //end loop on track collection
  } //end loop on the vector of track collections
  
  // no disablePreId_ switch, it is simpler to have an empty collection rather than no collection
  iEvent.put(std::move(output_preid),preidgsf_);
  if (produceCkfseed_)
    iEvent.put(std::move(output_nopre),preidckf_);
  if(producePreId_)
    {
      const edm::OrphanHandle<reco::PreIdCollection> preIdRefProd = iEvent.put(std::move(output_preidinfo),preidname_);
      // now make the Value Map, but only if one input collection
      if(tracksContainers_.size()==1)
    {
      Handle<TrackCollection> tkRefCollection ;
      iEvent.getByToken(tracksContainers_[0],tkRefCollection);
      fillPreIdRefValueMap(tkRefCollection,preIdRefProd,mapFiller);
      mapFiller.fill();
      iEvent.put(std::move(preIdMap_p),preidname_);
    }
    }
  // clear temporary maps
  refMap_.clear();
}

// intialize the cross-thread cache to hold gbr trees and resolution maps
namespace goodseedhelpers {
  HeavyObjectCache::HeavyObjectCache( const edm::ParameterSet& conf) {    
    // mvas
    const bool useTmva = conf.getUntrackedParameter<bool>("UseTMVA",false);
    
    if( useTmva ) {
      const std::string method_ = conf.getParameter<string>("TMVAMethod");
      std::array<edm::FileInPath,kMaxWeights> weights = {{ edm::FileInPath(conf.getParameter<string>("Weights1")),
                                                           edm::FileInPath(conf.getParameter<string>("Weights2")),
                                                           edm::FileInPath(conf.getParameter<string>("Weights3")),
                                                           edm::FileInPath(conf.getParameter<string>("Weights4")),
                                                           edm::FileInPath(conf.getParameter<string>("Weights5")),
                                                           edm::FileInPath(conf.getParameter<string>("Weights6")),
                                                           edm::FileInPath(conf.getParameter<string>("Weights7")),
                                                           edm::FileInPath(conf.getParameter<string>("Weights8")),
                                                           edm::FileInPath(conf.getParameter<string>("Weights9")) }};
            
      for(UInt_t j = 0; j < gbr.size(); ++j){
        TMVA::Reader reader("!Color:Silent");
                
        reader.AddVariable("NHits", &nhit);
        reader.AddVariable("NormChi", &chikfred);
        reader.AddVariable("dPtGSF", &dpt);
        reader.AddVariable("EoP", &eP);
        reader.AddVariable("ChiRatio", &chiRatio);
        reader.AddVariable("RedChi", &chired);
        reader.AddVariable("EcalDEta", &trk_ecalDeta);
        reader.AddVariable("EcalDPhi", &trk_ecalDphi);
        reader.AddVariable("pt", &pt);
        reader.AddVariable("eta", &eta);
        
        reader.BookMVA(method_, weights[j].fullPath().c_str());
        
        gbr[j].reset( new GBRForest( dynamic_cast<TMVA::MethodBDT*>( reader.FindMVA(method_) ) ) );
      }    
    }
  }
}

// ------------ method called once each job just before starting event loop  ------------
void 
GoodSeedProducer::beginRun(const edm::Run & run,
               const EventSetup& es)
{
  //Magnetic Field
  ESHandle<MagneticField> magneticField;
  es.get<IdealMagneticFieldRecord>().get(magneticField);
  B_=magneticField->inTesla(GlobalPoint(0,0,0));
  
  pfTransformer_.reset( new PFTrackTransformer(B_) );
  pfTransformer_->OnlyProp();
  
  //Resolution maps
    FileInPath ecalEtaMap(conf_.getParameter<string>("EtaMap"));
    FileInPath ecalPhiMap(conf_.getParameter<string>("PhiMap"));
    resMapEtaECAL_.reset( new PFResolutionMap("ECAL_eta",ecalEtaMap.fullPath().c_str()) );
    resMapPhiECAL_.reset( new PFResolutionMap("ECAL_phi",ecalPhiMap.fullPath().c_str()) );

  //read threshold
  FileInPath parFile(conf_.getParameter<string>("ThresholdFile"));
  ifstream ifs(parFile.fullPath().c_str());
  for (int iy=0;iy<81;++iy) ifs >> thr[iy];
}

int 
GoodSeedProducer::getBin(float eta, float pt){
  int ie=0;
  int ip=0;
  if (fabs(eta)<0.8) ie=0;
  else{ if (fabs(eta)<1.479) ie=1;
    else ie=2;
  }
  if (pt<6) ip=0;
  else {  if (pt<12) ip=1;     
    else ip=2;
  }
  int iep= ie*3+ip;
  LogDebug("GoodSeedProducer")<<"Track pt ="<<pt<<" eta="<<eta<<" bin="<<iep;
  return iep;
}

void GoodSeedProducer::fillPreIdRefValueMap( Handle<TrackCollection> tracks,
                         const edm::OrphanHandle<reco::PreIdCollection>& preidhandle,
                         edm::ValueMap<reco::PreIdRef>::Filler & filler)
{
  std::vector<reco::PreIdRef> values;

  unsigned ntracks=tracks->size();
  for(unsigned itrack=0;itrack<ntracks;++itrack)
   {
     reco::TrackRef theTrackRef(tracks,itrack);
     std::map<reco::TrackRef,unsigned>::const_iterator itcheck=refMap_.find(theTrackRef);
     if(itcheck==refMap_.end()) 
       {
     // the track has been early discarded
     values.push_back(reco::PreIdRef());
       }
     else
       {
     edm::Ref<reco::PreIdCollection> preIdRef(preidhandle,itcheck->second);
     values.push_back(preIdRef);
     //  std::cout << " Checking Refs " << (theTrackRef==preIdRef->trackRef()) << std::endl;
       }
   }
  filler.insert(tracks,values.begin(),values.end());
}

DEFINE_FWK_MODULE(GoodSeedProducer);