Reference Manual

Go to the documentation of this file.
00001 //
00002 // Package:         RecoTracker/FinalTrackSelectors
00003 // Class:           ConversionTrackMerger
00004 // 
00005 // Description:     Merger for ConversionTracks, adapted from SimpleTrackListMerger
00006 //
00007 // Original Author: J.Bendavid
00008 //
00009 // $Author: bendavid $
00010 // $Date: 2011/02/22 20:49:21 $
00011 // $Revision: 1.3 $
00012 //
00013 
00014 #include <memory>
00015 #include <string>
00016 #include <iostream>
00017 #include <cmath>
00018 #include <vector>
00019 
00020 #include "RecoEgamma/EgammaPhotonProducers/interface/ConversionTrackMerger.h"
00021 
00022 #include "DataFormats/TrackerRecHit2D/interface/SiStripMatchedRecHit2DCollection.h"
00023 #include "DataFormats/TrackerRecHit2D/interface/SiStripRecHit2DCollection.h"
00024 #include "DataFormats/TrackerRecHit2D/interface/SiStripRecHit1DCollection.h"
00025 #include "DataFormats/TrajectorySeed/interface/TrajectorySeedCollection.h"
00026 #include "DataFormats/TrackCandidate/interface/TrackCandidateCollection.h"
00027 
00028 #include "FWCore/Framework/interface/ESHandle.h"
00029 
00030 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00031 
00032 #include "Geometry/CommonDetUnit/interface/GeomDetUnit.h"
00033 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
00034 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
00035 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
00036 
00037 //#include "DataFormats/TrackReco/src/classes.h"
00038 
00039 #include "RecoTracker/TrackProducer/interface/ClusterRemovalRefSetter.h"
00040 #include "CommonTools/Statistics/interface/ChiSquaredProbability.h"
00041 
00042     
00043   ConversionTrackMerger::ConversionTrackMerger(edm::ParameterSet const& conf) : 
00044     conf_(conf)
00045   {
00046     produces<reco::ConversionTrackCollection>();
00047    
00048   }
00049 
00050 
00051   // Virtual destructor needed.
00052   ConversionTrackMerger::~ConversionTrackMerger() { }  
00053 
00054   // Functions that gets called by framework every event
00055   void ConversionTrackMerger::produce(edm::Event& e, const edm::EventSetup& es)
00056   {
00057     // retrieve producer name of input TrackCollection(s)
00058     std::string trackProducer1 = conf_.getParameter<std::string>("TrackProducer1");
00059     std::string trackProducer2 = conf_.getParameter<std::string>("TrackProducer2");
00060 
00061     double shareFrac =  conf_.getParameter<double>("ShareFrac");
00062     bool allowFirstHitShare = conf_.getParameter<bool>("allowFirstHitShare");
00063     bool checkCharge = conf_.getParameter<bool>("checkCharge");    
00064     double minProb = conf_.getParameter<double>("minProb");
00065     
00066     int outputPreferCollection = conf_.getParameter<int>("outputPreferCollection");
00067     int trackerOnlyPreferCollection = conf_.getParameter<int>("trackerOnlyPreferCollection");
00068     int arbitratedEcalSeededPreferCollection = conf_.getParameter<int>("arbitratedEcalSeededPreferCollection");    
00069     int arbitratedMergedPreferCollection = conf_.getParameter<int>("arbitratedMergedPreferCollection");
00070     int arbitratedMergedEcalGeneralPreferCollection = conf_.getParameter<int>("arbitratedMergedEcalGeneralPreferCollection");    
00071 
00072     // get Inputs 
00073     // if 1 input list doesn't exist, make an empty list, issue a warning, and continue
00074     // this allows ConversionTrackMerger to be used as a cleaner only if handed just one list
00075     // if both input lists don't exist, will issue 2 warnings and generate an empty output collection
00076     //
00077     const reco::ConversionTrackCollection *TC1 = 0;
00078     static const reco::ConversionTrackCollection s_empty1, s_empty2;
00079     edm::Handle<reco::ConversionTrackCollection> trackCollection1;
00080     e.getByLabel(trackProducer1, trackCollection1);
00081     if (trackCollection1.isValid()) {
00082       TC1 = trackCollection1.product();
00083       //std::cout << "1st collection " << trackProducer1 << " has "<< TC1->size() << " tracks" << std::endl ;
00084     } else {
00085       TC1 = &s_empty1;
00086       edm::LogWarning("ConversionTrackMerger") << "1st TrackCollection " << trackProducer1 << " not found; will only clean 2nd TrackCollection " << trackProducer2 ;
00087     }
00088     reco::ConversionTrackCollection tC1 = *TC1;
00089 
00090     const reco::ConversionTrackCollection *TC2 = 0;
00091     edm::Handle<reco::ConversionTrackCollection> trackCollection2;
00092     e.getByLabel(trackProducer2, trackCollection2);
00093     if (trackCollection2.isValid()) {
00094       TC2 = trackCollection2.product();
00095       //std::cout << "2nd collection " << trackProducer2 << " has "<< TC2->size() << " tracks" << std::endl ;
00096     } else {
00097         TC2 = &s_empty2;
00098         edm::LogWarning("ConversionTrackMerger") << "2nd TrackCollection " << trackProducer2 << " not found; will only clean 1st TrackCollection " << trackProducer1 ;
00099     }
00100     reco::ConversionTrackCollection tC2 = *TC2;
00101 
00102     // Step B: create empty output collection
00103     outputTrks = std::auto_ptr<reco::ConversionTrackCollection>(new reco::ConversionTrackCollection);
00104     int i;
00105 
00106     std::vector<int> selected1; for (unsigned int i=0; i<tC1.size(); ++i){selected1.push_back(1);}
00107     std::vector<int> selected2; for (unsigned int i=0; i<tC2.size(); ++i){selected2.push_back(1);}
00108 
00109    
00110    std::map<reco::ConversionTrackCollection::const_iterator, std::vector<const TrackingRecHit*> > rh1;
00111    std::map<reco::ConversionTrackCollection::const_iterator, std::vector<const TrackingRecHit*> > rh2;
00112    for (reco::ConversionTrackCollection::const_iterator track=tC1.begin(); track!=tC1.end(); ++track){
00113      trackingRecHit_iterator itB = track->track()->recHitsBegin();
00114      trackingRecHit_iterator itE = track->track()->recHitsEnd();
00115      for (trackingRecHit_iterator it = itB;  it != itE; ++it) { 
00116        const TrackingRecHit* hit = &(**it);
00117        rh1[track].push_back(hit);
00118      }
00119    }
00120    for (reco::ConversionTrackCollection::const_iterator track=tC2.begin(); track!=tC2.end(); ++track){
00121      trackingRecHit_iterator jtB = track->track()->recHitsBegin();
00122      trackingRecHit_iterator jtE = track->track()->recHitsEnd();
00123      for (trackingRecHit_iterator jt = jtB;  jt != jtE; ++jt) { 
00124        const TrackingRecHit* hit = &(**jt);
00125        rh2[track].push_back(hit);
00126      }
00127    }
00128 
00129    if ( (0<tC1.size())&&(0<tC2.size()) ){
00130     i=-1;
00131     for (reco::ConversionTrackCollection::iterator track=tC1.begin(); track!=tC1.end(); ++track){
00132       i++; 
00133         
00134       //clear flags if preferCollection was set to 0
00135       selected1[i] = selected1[i] && outputPreferCollection!=0;
00136       track->setIsTrackerOnly ( track->isTrackerOnly() &&  trackerOnlyPreferCollection!=0 );
00137       track->setIsArbitratedEcalSeeded( track->isArbitratedEcalSeeded() &&  arbitratedEcalSeededPreferCollection!=0 );
00138       track->setIsArbitratedMerged( track->isArbitratedMerged() && arbitratedMergedPreferCollection!=0 );
00139       track->setIsArbitratedMergedEcalGeneral( track->isArbitratedMergedEcalGeneral() && arbitratedMergedEcalGeneralPreferCollection!=0 );
00140       
00141 
00142       std::vector<const TrackingRecHit*>& iHits = rh1[track]; 
00143       unsigned nh1 = iHits.size();
00144       int j=-1;
00145       for (reco::ConversionTrackCollection::iterator track2=tC2.begin(); track2!=tC2.end(); ++track2){
00146         j++;
00147 
00148         //clear flags if preferCollection was set to 0
00149         selected2[j] = selected2[j] && outputPreferCollection!=0;
00150         track2->setIsTrackerOnly ( track2->isTrackerOnly() &&  trackerOnlyPreferCollection!=0 );
00151         track2->setIsArbitratedEcalSeeded( track2->isArbitratedEcalSeeded() &&  arbitratedEcalSeededPreferCollection!=0 );
00152         track2->setIsArbitratedMerged( track2->isArbitratedMerged() && arbitratedMergedPreferCollection!=0 );
00153         track2->setIsArbitratedMergedEcalGeneral( track2->isArbitratedMergedEcalGeneral() && arbitratedMergedEcalGeneralPreferCollection!=0 );
00154 
00155         std::vector<const TrackingRecHit*>& jHits = rh2[track2]; 
00156         unsigned nh2 = jHits.size();
00157         int noverlap=0;
00158         int firstoverlap=0;
00159         for ( unsigned ih=0; ih<nh1; ++ih ) { 
00160           const TrackingRecHit* it = iHits[ih];
00161           if (it->isValid()){
00162             int jj=-1;
00163             for ( unsigned jh=0; jh<nh2; ++jh ) { 
00164               const TrackingRecHit* jt = jHits[jh];
00165               jj++;
00166               if (jt->isValid()){           
00167                 if ( it->sharesInput(jt,TrackingRecHit::some) ) {
00168                   noverlap++;
00169                   if ( allowFirstHitShare && ( ih == 0 ) && ( jh == 0 ) ) firstoverlap=1;
00170                 }
00171               }
00172             }
00173           }
00174         }
00175         int nhit1 = track->track()->numberOfValidHits();
00176         int nhit2 = track2->track()->numberOfValidHits();
00177         //if (noverlap>0) printf("noverlap = %i, firstoverlap = %i, nhit1 = %i, nhit2 = %i, algo1 = %i, algo2 = %i, q1 = %i, q2 = %i\n",noverlap,firstoverlap,nhit1,nhit2,track->track()->algo(),track2->track()->algo(),track->track()->charge(),track2->track()->charge());
00178         //std::cout << " trk1 trk2 nhits1 nhits2 nover " << i << " " << j << " " << track->numberOfValidHits() << " "  << track2->numberOfValidHits() << " " << noverlap << " " << fi << " " << fj  <<std::endl;
00179         if ( (noverlap-firstoverlap) > (std::min(nhit1,nhit2)-firstoverlap)*shareFrac && (!checkCharge || track->track()->charge()*track2->track()->charge()>0) ) {
00180           //printf("overlapping tracks\n");
00181          //printf ("ndof1 = %5f, chisq1 = %5f, ndof2 = %5f, chisq2 = %5f\n",track->track()->ndof(),track->track()->chi2(),track2->track()->ndof(),track2->track()->chi2());
00182           
00183           double probFirst = ChiSquaredProbability(track->track()->chi2(), track->track()->ndof());
00184           double probSecond = ChiSquaredProbability(track2->track()->chi2(), track2->track()->ndof());
00185 
00186           //arbitrate by number of hits and reduced chisq
00187           bool keepFirst = ( nhit1>nhit2 || (nhit1==nhit2 && track->track()->normalizedChi2()<track2->track()->normalizedChi2()) );
00188 
00189           //override decision in case one track is radically worse quality than the other
00190           keepFirst |= (probFirst>minProb && probSecond<=minProb);
00191           keepFirst &= !(probFirst<=minProb && probSecond>minProb);
00192 
00193           bool keepSecond = !keepFirst;
00194                     
00195           //set flags based on arbitration decision and precedence settings
00196 
00197           selected1[i] =            selected1[i]            && ( (keepFirst && outputPreferCollection==3) || outputPreferCollection==-1 || outputPreferCollection==1 );
00198           track->setIsTrackerOnly ( track->isTrackerOnly() && ( (keepFirst && trackerOnlyPreferCollection==3) || trackerOnlyPreferCollection==-1 || trackerOnlyPreferCollection==1 ) );
00199           track->setIsArbitratedEcalSeeded( track->isArbitratedEcalSeeded() &&  ( (keepFirst && arbitratedEcalSeededPreferCollection==3) || arbitratedEcalSeededPreferCollection==-1 || arbitratedEcalSeededPreferCollection==1 ) );
00200           track->setIsArbitratedMerged( track->isArbitratedMerged() && ( (keepFirst && arbitratedMergedPreferCollection==3) || arbitratedMergedPreferCollection==-1 || arbitratedMergedPreferCollection==1 ) );
00201           track->setIsArbitratedMergedEcalGeneral( track->isArbitratedMergedEcalGeneral() && ( (keepFirst && arbitratedMergedEcalGeneralPreferCollection==3) || arbitratedMergedEcalGeneralPreferCollection==-1 || arbitratedMergedEcalGeneralPreferCollection==1 ) );
00202           
00203           selected2[j] =             selected2[j]            && ( (keepSecond && outputPreferCollection==3) || outputPreferCollection==-1 || outputPreferCollection==2 );
00204           track2->setIsTrackerOnly ( track2->isTrackerOnly() && ( (keepSecond && trackerOnlyPreferCollection==3) || trackerOnlyPreferCollection==-1 || trackerOnlyPreferCollection==2 ) );
00205           track2->setIsArbitratedEcalSeeded( track2->isArbitratedEcalSeeded() &&  ( (keepSecond && arbitratedEcalSeededPreferCollection==3) || arbitratedEcalSeededPreferCollection==-1 || arbitratedEcalSeededPreferCollection==2 ) );
00206           track2->setIsArbitratedMerged( track2->isArbitratedMerged() && ( (keepSecond && arbitratedMergedPreferCollection==3) || arbitratedMergedPreferCollection==-1 || arbitratedMergedPreferCollection==2 ) );
00207           track2->setIsArbitratedMergedEcalGeneral( track2->isArbitratedMergedEcalGeneral() && ( (keepSecond && arbitratedMergedEcalGeneralPreferCollection==3) || arbitratedMergedEcalGeneralPreferCollection==-1 || arbitratedMergedEcalGeneralPreferCollection==2 ) );
00208           
00209         }//end got a duplicate
00210       }//end track2 loop
00211     }//end track loop
00212    }//end more than 1 track
00213 
00214   //
00215   //  output selected tracks - if any
00216   //
00217    
00218    if ( 0<tC1.size() ){
00219      i=0;
00220      for (reco::ConversionTrackCollection::const_iterator track=tC1.begin(); track!=tC1.end(); 
00221           ++track, ++i){
00222       //don't store tracks rejected as duplicates
00223       if (!selected1[i]){
00224         continue;
00225       }
00226       //fill the TrackCollection
00227       outputTrks->push_back(*track);      
00228     }//end faux loop over tracks
00229    }//end more than 0 track
00230 
00231    //Fill the trajectories, etc. for 1st collection
00232  
00233    
00234    if ( 0<tC2.size() ){
00235     i=0;
00236     for (reco::ConversionTrackCollection::const_iterator track=tC2.begin(); track!=tC2.end();
00237          ++track, ++i){
00238       //don't store tracks rejected as duplicates
00239       if (!selected2[i]){
00240         continue;
00241       }
00242       //fill the TrackCollection
00243       outputTrks->push_back(*track);
00244     }//end faux loop over tracks
00245    }//end more than 0 track
00246   
00247     e.put(outputTrks);
00248     return;
00249 
00250   }//end produce
/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_6_2_5/src/RecoEgamma/EgammaPhotonProducers/src/ConversionTrackMerger.cc
