Reference Manual

Go to the documentation of this file.
00001 // -*- C++ -*-
00002 //
00003 // Package:    PFTracking
00004 // Class:      PFElecTkProducer
00005 // 
00006 // Original Author:  Michele Pioppi
00007 //         Created:  Tue Jan 23 15:26:39 CET 2007
00008 
00009 
00010 
00011 // system include files
00012 #include <memory>
00013 
00014 // user include files
00015 #include "RecoParticleFlow/PFTracking/interface/PFElecTkProducer.h"
00016 #include "RecoParticleFlow/PFTracking/interface/PFTrackTransformer.h"
00017 #include "RecoParticleFlow/PFTracking/interface/ConvBremPFTrackFinder.h"
00018 #include "DataFormats/GsfTrackReco/interface/GsfTrackFwd.h"
00019 #include "MagneticField/Engine/interface/MagneticField.h"
00020 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
00021 #include "FWCore/Framework/interface/ESHandle.h"
00022 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00023 #include "TrackingTools/PatternTools/interface/Trajectory.h"
00024 #include "DataFormats/TrajectorySeed/interface/TrajectorySeed.h"
00025 #include "DataFormats/EgammaCandidates/interface/GsfElectronFwd.h"
00026 #include "DataFormats/EgammaCandidates/interface/GsfElectron.h"
00027 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
00028 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
00029 #include "MagneticField/Engine/interface/MagneticField.h"
00030 #include "DataFormats/ParticleFlowReco/interface/PFClusterFwd.h"
00031 #include "DataFormats/ParticleFlowReco/interface/PFCluster.h"
00032 #include "TrackingTools/TransientTrack/interface/TransientTrackBuilder.h"
00033 #include "TrackingTools/Records/interface/TransientTrackRecord.h"
00034 #include "DataFormats/VertexReco/interface/Vertex.h"
00035 #include "DataFormats/VertexReco/interface/VertexFwd.h"
00036 #include "DataFormats/ParticleFlowReco/interface/PFDisplacedTrackerVertex.h"
00037 #include "DataFormats/ParticleFlowReco/interface/PFDisplacedVertexFwd.h"
00038 #include "DataFormats/ParticleFlowReco/interface/PFDisplacedVertex.h"
00039 #include "DataFormats/ParticleFlowReco/interface/PFConversionFwd.h"
00040 #include "DataFormats/ParticleFlowReco/interface/PFConversion.h"
00041 #include "DataFormats/ParticleFlowReco/interface/PFV0Fwd.h"
00042 #include "DataFormats/ParticleFlowReco/interface/PFV0.h"
00043 #include "RecoParticleFlow/PFClusterTools/interface/LinkByRecHit.h"
00044 #include "RecoParticleFlow/PFClusterTools/interface/ClusterClusterMapping.h"
00045 
00046 #include "TMath.h"
00047 using namespace std;
00048 using namespace edm;
00049 using namespace reco;
00050 PFElecTkProducer::PFElecTkProducer(const ParameterSet& iConfig):
00051   conf_(iConfig),
00052   pfTransformer_(0),
00053   convBremFinder_(0)
00054 {
00055   LogInfo("PFElecTkProducer")<<"PFElecTkProducer started";
00056 
00057   gsfTrackLabel_ = iConfig.getParameter<InputTag>
00058     ("GsfTrackModuleLabel");
00059 
00060   pfTrackLabel_ = iConfig.getParameter<InputTag>
00061     ("PFRecTrackLabel");
00062 
00063   primVtxLabel_ = iConfig.getParameter<InputTag>
00064     ("PrimaryVertexLabel");
00065 
00066   pfEcalClusters_ = iConfig.getParameter<InputTag>
00067     ("PFEcalClusters");
00068 
00069   pfNuclear_ = iConfig.getParameter<InputTag>
00070     ("PFNuclear");
00071   
00072   pfConv_ = iConfig.getParameter<InputTag>
00073     ("PFConversions");
00074   
00075   pfV0_ = iConfig.getParameter<InputTag>
00076     ("PFV0");
00077 
00078   useNuclear_ = iConfig.getParameter<bool>("useNuclear");
00079   useConversions_ = iConfig.getParameter<bool>("useConversions");
00080   useV0_ = iConfig.getParameter<bool>("useV0");
00081   debugGsfCleaning_ = iConfig.getParameter<bool>("debugGsfCleaning");
00082 
00083   produces<GsfPFRecTrackCollection>();
00084   produces<GsfPFRecTrackCollection>( "Secondary" ).setBranchAlias( "secondary" );
00085 
00086 
00087   trajinev_ = iConfig.getParameter<bool>("TrajInEvents");
00088   modemomentum_ = iConfig.getParameter<bool>("ModeMomentum");
00089   applySel_ = iConfig.getParameter<bool>("applyEGSelection");
00090   applyGsfClean_ = iConfig.getParameter<bool>("applyGsfTrackCleaning");
00091   applyAngularGsfClean_ = iConfig.getParameter<bool>("applyAlsoGsfAngularCleaning");
00092   detaCutGsfClean_ =  iConfig.getParameter<double>("maxDEtaGsfAngularCleaning");
00093   dphiCutGsfClean_ =  iConfig.getParameter<double>("maxDPhiBremTangGsfAngularCleaning");
00094   useFifthStepForTrackDriven_ = iConfig.getParameter<bool>("useFifthStepForTrackerDrivenGsf");
00095   useFifthStepForEcalDriven_ = iConfig.getParameter<bool>("useFifthStepForEcalDrivenGsf");
00096   maxPtConvReco_ = iConfig.getParameter<double>("MaxConvBremRecoPT");
00097   detaGsfSC_ = iConfig.getParameter<double>("MinDEtaGsfSC");
00098   dphiGsfSC_ = iConfig.getParameter<double>("MinDPhiGsfSC");
00099   SCEne_ = iConfig.getParameter<double>("MinSCEnergy");
00100   
00101   // set parameter for convBremFinder
00102   useConvBremFinder_ =     iConfig.getParameter<bool>("useConvBremFinder");
00103   mvaConvBremFinderID_
00104     = iConfig.getParameter<double>("pf_convBremFinderID_mvaCut");
00105   
00106   string mvaWeightFileConvBrem
00107     = iConfig.getParameter<string>("pf_convBremFinderID_mvaWeightFile");
00108   
00109   
00110   if(useConvBremFinder_) 
00111     path_mvaWeightFileConvBrem_ = edm::FileInPath ( mvaWeightFileConvBrem.c_str() ).fullPath();
00112 
00113 }
00114 
00115 
00116 PFElecTkProducer::~PFElecTkProducer()
00117 {
00118  
00119   delete pfTransformer_;
00120   delete convBremFinder_;
00121 }
00122 
00123 
00124 //
00125 // member functions
00126 //
00127 
00128 // ------------ method called to produce the data  ------------
00129 void
00130 PFElecTkProducer::produce(Event& iEvent, const EventSetup& iSetup)
00131 {
00132   LogDebug("PFElecTkProducer")<<"START event: "<<iEvent.id().event()
00133                               <<" in run "<<iEvent.id().run();
00134 
00135   //create the empty collections 
00136   auto_ptr< GsfPFRecTrackCollection > 
00137     gsfPFRecTrackCollection(new GsfPFRecTrackCollection);
00138 
00139   
00140   auto_ptr< GsfPFRecTrackCollection > 
00141     gsfPFRecTrackCollectionSecondary(new GsfPFRecTrackCollection);
00142 
00143   //read collections of tracks
00144   Handle<GsfTrackCollection> gsftrackscoll;
00145   iEvent.getByLabel(gsfTrackLabel_,gsftrackscoll);
00146 
00147   //read collections of trajectories
00148   Handle<vector<Trajectory> > TrajectoryCollection;
00149  
00150   //read pfrectrack collection
00151   Handle<PFRecTrackCollection> thePfRecTrackCollection;
00152   iEvent.getByLabel(pfTrackLabel_,thePfRecTrackCollection);
00153   const PFRecTrackCollection& PfRTkColl = *(thePfRecTrackCollection.product());
00154 
00155   // PFClusters
00156   Handle<PFClusterCollection> theECPfClustCollection;
00157   iEvent.getByLabel(pfEcalClusters_,theECPfClustCollection);
00158   const PFClusterCollection& theEcalClusters = *(theECPfClustCollection.product());
00159 
00160   //Primary Vertexes
00161   Handle<reco::VertexCollection> thePrimaryVertexColl;
00162   iEvent.getByLabel(primVtxLabel_,thePrimaryVertexColl);
00163 
00164 
00165 
00166   // Displaced Vertex
00167   Handle< reco::PFDisplacedTrackerVertexCollection > pfNuclears; 
00168   if( useNuclear_ ) {
00169     bool found = iEvent.getByLabel(pfNuclear_, pfNuclears);
00170     
00171     
00172     if(!found )
00173       LogError("PFElecTkProducer")<<" cannot get PFNuclear : "
00174                                   <<  pfNuclear_
00175                                   << " please set useNuclear=False in RecoParticleFlow/PFTracking/python/pfTrackElec_cfi.py" << endl;
00176   }
00177 
00178   // Conversions 
00179   Handle< reco::PFConversionCollection > pfConversions;
00180   if( useConversions_ ) {
00181     bool found = iEvent.getByLabel(pfConv_,pfConversions);
00182     if(!found )
00183       LogError("PFElecTkProducer")<<" cannot get PFConversions : "
00184                                   << pfConv_ 
00185                                   << " please set useConversions=False in RecoParticleFlow/PFTracking/python/pfTrackElec_cfi.py" << endl;
00186   }
00187 
00188   // V0
00189   Handle< reco::PFV0Collection > pfV0;
00190   if( useV0_ ) {
00191     bool found = iEvent.getByLabel(pfV0_, pfV0);
00192     
00193     if(!found )
00194       LogError("PFElecTkProducer")<<" cannot get PFV0 : "
00195                                   << pfV0_
00196                                   << " please set useV0=False  RecoParticleFlow/PFTracking/python/pfTrackElec_cfi.py" << endl;
00197   }
00198   
00199   
00200 
00201   GsfTrackCollection gsftracks = *(gsftrackscoll.product());    
00202   vector<Trajectory> tjvec(0);
00203   if (trajinev_){
00204     bool foundTraj = iEvent.getByLabel(gsfTrackLabel_,TrajectoryCollection); 
00205     if(!foundTraj) 
00206       LogError("PFElecTkProducer")
00207         <<" cannot get Trajectories of : "
00208         <<  gsfTrackLabel_
00209         << " please set TrajInEvents = False in RecoParticleFlow/PFTracking/python/pfTrackElec_cfi.py" << endl;
00210     
00211     tjvec= *(TrajectoryCollection.product());
00212   }
00213   
00214   
00215   vector<reco::GsfPFRecTrack> selGsfPFRecTracks;
00216   vector<reco::GsfPFRecTrack> primaryGsfPFRecTracks;
00217   std::map<unsigned int, std::vector<reco::GsfPFRecTrack> >  GsfPFMap;
00218   
00219 
00220   for (unsigned int igsf=0; igsf<gsftracks.size();igsf++) {
00221     
00222     GsfTrackRef trackRef(gsftrackscoll, igsf);
00223 
00224     int kf_ind=FindPfRef(PfRTkColl,gsftracks[igsf],false);
00225     
00226     if (kf_ind>=0) {
00227       
00228       PFRecTrackRef kf_ref(thePfRecTrackCollection,
00229                            kf_ind);
00230 
00231       // remove fifth step tracks
00232       if( useFifthStepForEcalDriven_ == false
00233           || useFifthStepForTrackDriven_ == false) {
00234         bool isFifthStepTrack = isFifthStep(kf_ref);
00235         bool isEcalDriven = true;
00236         bool isTrackerDriven = true;
00237         
00238         if (&(*trackRef->seedRef())==0) {
00239           isEcalDriven = false;
00240           isTrackerDriven = false;
00241         }
00242         else {
00243           ElectronSeedRef SeedRef= trackRef->extra()->seedRef().castTo<ElectronSeedRef>();
00244           if(SeedRef->caloCluster().isNull())
00245             isEcalDriven = false;
00246           if(SeedRef->ctfTrack().isNull())
00247             isTrackerDriven = false;
00248         }
00249         //note: the same track could be both ecalDriven and trackerDriven
00250         if(isFifthStepTrack && 
00251            isEcalDriven &&
00252            isTrackerDriven == false &&
00253            useFifthStepForEcalDriven_ == false) {
00254           continue;
00255         }
00256 
00257         if(isFifthStepTrack && 
00258            isTrackerDriven  && 
00259            isEcalDriven == false && 
00260            useFifthStepForTrackDriven_ == false) {
00261           continue;
00262         }
00263 
00264         if(isFifthStepTrack && 
00265            isTrackerDriven && 
00266            isEcalDriven && 
00267            useFifthStepForTrackDriven_ == false &&
00268            useFifthStepForEcalDriven_ == false) {
00269           continue;
00270         }
00271       }
00272       
00273       pftrack_=GsfPFRecTrack( gsftracks[igsf].charge(), 
00274                               reco::PFRecTrack::GSF, 
00275                               igsf, trackRef,
00276                               kf_ref);
00277     } else  {
00278       PFRecTrackRef dummyRef;
00279       pftrack_=GsfPFRecTrack( gsftracks[igsf].charge(), 
00280                               reco::PFRecTrack::GSF, 
00281                               igsf, trackRef,
00282                               dummyRef);
00283     }
00284     
00285     
00286     bool validgsfbrem = false;
00287     if(trajinev_) {
00288       validgsfbrem = pfTransformer_->addPointsAndBrems(pftrack_, 
00289                                                        gsftracks[igsf], 
00290                                                        tjvec[igsf],
00291                                                        modemomentum_);
00292     } else {
00293       validgsfbrem = pfTransformer_->addPointsAndBrems(pftrack_, 
00294                                                        gsftracks[igsf], 
00295                                                        mtsTransform_);
00296     }
00297     
00298     bool passSel = true;
00299     if(applySel_) 
00300       passSel = applySelection(gsftracks[igsf]);      
00301 
00302     if(validgsfbrem && passSel) 
00303       selGsfPFRecTracks.push_back(pftrack_);
00304   }
00305   
00306 
00307   unsigned int count_primary = 0;
00308   if(selGsfPFRecTracks.size() > 0) {
00309     for(unsigned int ipfgsf=0; ipfgsf<selGsfPFRecTracks.size();ipfgsf++) {
00310       
00311       vector<unsigned int> secondaries(0);
00312       secondaries.clear();
00313       bool keepGsf = true;
00314 
00315       if(applyGsfClean_) {
00316         keepGsf = resolveGsfTracks(selGsfPFRecTracks,ipfgsf,secondaries,theEcalClusters);
00317       }
00318  
00319       //is primary? 
00320       if(keepGsf == true) {
00321 
00322         // Find kf tracks from converted brem photons
00323         if(convBremFinder_->foundConvBremPFRecTrack(thePfRecTrackCollection,thePrimaryVertexColl,
00324                                                     pfNuclears,pfConversions,pfV0,
00325                                                     useNuclear_,useConversions_,useV0_,
00326                                                     theEcalClusters,selGsfPFRecTracks[ipfgsf])) {
00327           const vector<PFRecTrackRef>& convBremPFRecTracks(convBremFinder_->getConvBremPFRecTracks());
00328           for(unsigned int ii = 0; ii<convBremPFRecTracks.size(); ii++) {
00329             selGsfPFRecTracks[ipfgsf].addConvBremPFRecTrackRef(convBremPFRecTracks[ii]);
00330           }
00331         }
00332         
00333         // save primaries gsf tracks
00334         //      gsfPFRecTrackCollection->push_back(selGsfPFRecTracks[ipfgsf]);
00335         primaryGsfPFRecTracks.push_back(selGsfPFRecTracks[ipfgsf]);
00336         
00337         
00338         // NOTE:: THE TRACKID IS USED TO LINK THE PRIMARY GSF TRACK. THIS NEEDS 
00339         // TO BE CHANGED AS SOON AS IT IS POSSIBLE TO CHANGE DATAFORMATS
00340         // A MODIFICATION HERE IMPLIES A MODIFICATION IN PFBLOCKALGO.CC/H
00341         unsigned int primGsfIndex = selGsfPFRecTracks[ipfgsf].trackId();
00342         vector<reco::GsfPFRecTrack> trueGsfPFRecTracks;
00343         if(secondaries.size() > 0) {
00344           // loop on secondaries gsf tracks (from converted brems)
00345           for(unsigned int isecpfgsf=0; isecpfgsf<secondaries.size();isecpfgsf++) {
00346             
00347             PFRecTrackRef refsecKF =  selGsfPFRecTracks[(secondaries[isecpfgsf])].kfPFRecTrackRef();
00348             
00349             unsigned int secGsfIndex = selGsfPFRecTracks[(secondaries[isecpfgsf])].trackId();
00350             GsfTrackRef secGsfRef = selGsfPFRecTracks[(secondaries[isecpfgsf])].gsfTrackRef();
00351 
00352             if(refsecKF.isNonnull()) {
00353               // NOTE::IT SAVED THE TRACKID OF THE PRIMARY!!! THIS IS USED IN PFBLOCKALGO.CC/H
00354               secpftrack_= GsfPFRecTrack( gsftracks[secGsfIndex].charge(), 
00355                                           reco::PFRecTrack::GSF, 
00356                                           primGsfIndex, secGsfRef,
00357                                           refsecKF);
00358             }
00359             else{
00360               PFRecTrackRef dummyRef;
00361               // NOTE::IT SAVED THE TRACKID OF THE PRIMARY!!! THIS IS USED IN PFBLOCKALGO.CC/H
00362               secpftrack_= GsfPFRecTrack( gsftracks[secGsfIndex].charge(), 
00363                                           reco::PFRecTrack::GSF, 
00364                                           primGsfIndex, secGsfRef,
00365                                           dummyRef);
00366             }
00367 
00368             bool validgsfbrem = false;
00369             if(trajinev_) {
00370               validgsfbrem = pfTransformer_->addPointsAndBrems(secpftrack_,
00371                                                                gsftracks[secGsfIndex], 
00372                                                                tjvec[secGsfIndex],
00373                                                                modemomentum_);
00374             } else {
00375               validgsfbrem = pfTransformer_->addPointsAndBrems(secpftrack_,
00376                                                                gsftracks[secGsfIndex], 
00377                                                                mtsTransform_);
00378             }         
00379 
00380             if(validgsfbrem) {
00381               gsfPFRecTrackCollectionSecondary->push_back(secpftrack_);
00382               trueGsfPFRecTracks.push_back(secpftrack_);
00383             }       
00384           }
00385         }
00386         GsfPFMap.insert(pair<unsigned int,std::vector<reco::GsfPFRecTrack> >(count_primary,trueGsfPFRecTracks));
00387         trueGsfPFRecTracks.clear();
00388         count_primary++;
00389       }
00390     }
00391   }
00392   
00393   
00394   const edm::OrphanHandle<GsfPFRecTrackCollection> gsfPfRefProd = 
00395     iEvent.put(gsfPFRecTrackCollectionSecondary,"Secondary");
00396   
00397   
00398   //now the secondary GsfPFRecTracks are in the event, the Ref can be created
00399   createGsfPFRecTrackRef(gsfPfRefProd,primaryGsfPFRecTracks,GsfPFMap);
00400   
00401   for(unsigned int iGSF = 0; iGSF<primaryGsfPFRecTracks.size();iGSF++){
00402     gsfPFRecTrackCollection->push_back(primaryGsfPFRecTracks[iGSF]);
00403   }
00404   iEvent.put(gsfPFRecTrackCollection);
00405 
00406   selGsfPFRecTracks.clear();
00407   GsfPFMap.clear();
00408   primaryGsfPFRecTracks.clear();
00409 }
00410 
00411 // create the secondary GsfPFRecTracks Ref
00412 void
00413 PFElecTkProducer::createGsfPFRecTrackRef(const edm::OrphanHandle<reco::GsfPFRecTrackCollection>& gsfPfHandle,
00414                                          std::vector<reco::GsfPFRecTrack>& gsfPFRecTrackPrimary,
00415                                          const std::map<unsigned int, std::vector<reco::GsfPFRecTrack> >& MapPrimSec) {
00416   unsigned int cgsf=0;
00417   unsigned int csecgsf=0;
00418   for (std::map<unsigned int, std::vector<reco::GsfPFRecTrack> >::const_iterator igsf = MapPrimSec.begin();
00419        igsf != MapPrimSec.end(); igsf++,cgsf++) {
00420     vector<reco::GsfPFRecTrack> SecGsfPF = igsf->second;
00421     for (unsigned int iSecGsf=0; iSecGsf < SecGsfPF.size(); iSecGsf++) {
00422       edm::Ref<reco::GsfPFRecTrackCollection> refgprt(gsfPfHandle,csecgsf);
00423       gsfPFRecTrackPrimary[cgsf].addConvBremGsfPFRecTrackRef(refgprt);
00424       ++csecgsf;
00425     }
00426   }
00427 
00428   return;
00429 }
00430 // ------------- method for find the corresponding kf pfrectrack ---------------------
00431 int
00432 PFElecTkProducer::FindPfRef(const reco::PFRecTrackCollection  & PfRTkColl, 
00433                             reco::GsfTrack gsftk,
00434                             bool otherColl){
00435 
00436 
00437   if (&(*gsftk.seedRef())==0) return -1;
00438   ElectronSeedRef ElSeedRef=gsftk.extra()->seedRef().castTo<ElectronSeedRef>();
00439   //CASE 1 ELECTRONSEED DOES NOT HAVE A REF TO THE CKFTRACK
00440   if (ElSeedRef->ctfTrack().isNull()){
00441     reco::PFRecTrackCollection::const_iterator pft=PfRTkColl.begin();
00442     reco::PFRecTrackCollection::const_iterator pftend=PfRTkColl.end();
00443     unsigned int i_pf=0;
00444     int ibest=-1;
00445     unsigned int ish_max=0;
00446     float dr_min=1000;
00447     //SEARCH THE PFRECTRACK THAT SHARES HITS WITH THE ELECTRON SEED
00448     // Here the cpu time can be improved. 
00449     for(;pft!=pftend;++pft){
00450       unsigned int ish=0;
00451       
00452       float dph= fabs(pft->trackRef()->phi()-gsftk.phi()); 
00453       if (dph>TMath::Pi()) dph-= TMath::TwoPi();
00454       float det=fabs(pft->trackRef()->eta()-gsftk.eta());
00455       float dr =sqrt(dph*dph+det*det);  
00456       
00457       trackingRecHit_iterator  hhit=
00458         pft->trackRef()->recHitsBegin();
00459       trackingRecHit_iterator  hhit_end=
00460         pft->trackRef()->recHitsEnd();
00461       
00462     
00463       
00464       for(;hhit!=hhit_end;++hhit){
00465         if (!(*hhit)->isValid()) continue;
00466         TrajectorySeed::const_iterator hit=
00467           gsftk.seedRef()->recHits().first;
00468         TrajectorySeed::const_iterator hit_end=
00469           gsftk.seedRef()->recHits().second;
00470         for(;hit!=hit_end;++hit){
00471           if (!(hit->isValid())) continue;
00472           if((*hhit)->sharesInput(&*(hit),TrackingRecHit::all))  ish++; 
00473         //   if((hit->geographicalId()==(*hhit)->geographicalId())&&
00474         //     (((*hhit)->localPosition()-hit->localPosition()).mag()<0.01)) ish++;
00475         }       
00476         
00477       }
00478       
00479 
00480       if ((ish>ish_max)||
00481           ((ish==ish_max)&&(dr<dr_min))){
00482         ish_max=ish;
00483         dr_min=dr;
00484         ibest=i_pf;
00485       }
00486       
00487    
00488     
00489       i_pf++;
00490     }
00491     if (ibest<0) return -1;
00492     
00493     if((ish_max==0) || (dr_min>0.05))return -1;
00494     if(otherColl && (ish_max==0)) return -1;
00495     return ibest;
00496   }
00497   else{
00498     //ELECTRON SEED HAS A REFERENCE
00499    
00500     reco::PFRecTrackCollection::const_iterator pft=PfRTkColl.begin();
00501     reco::PFRecTrackCollection::const_iterator pftend=PfRTkColl.end();
00502     unsigned int i_pf=0;
00503     
00504     for(;pft!=pftend;++pft){
00505       //REF COMPARISON
00506       if (pft->trackRef()==ElSeedRef->ctfTrack()){
00507         return i_pf;
00508       }
00509       i_pf++;
00510     }
00511   }
00512   return -1;
00513 }
00514 bool PFElecTkProducer::isFifthStep(reco::PFRecTrackRef pfKfTrack) {
00515 
00516   bool isFithStep = false;
00517   
00518 
00519   TrackRef kfref = pfKfTrack->trackRef();
00520   unsigned int Algo = 0; 
00521   switch (kfref->algo()) {
00522   case TrackBase::undefAlgorithm:
00523   case TrackBase::ctf:
00524   case TrackBase::iter0:
00525   case TrackBase::iter1:
00526   case TrackBase::iter2:
00527     Algo = 0;
00528     break;
00529   case TrackBase::iter3:
00530     Algo = 1;
00531     break;
00532   case TrackBase::iter4:
00533     Algo = 2;
00534     break;
00535   case TrackBase::iter5:
00536     Algo = 3;
00537     break;
00538   case TrackBase::iter6:
00539     Algo = 4;
00540     break;
00541   default:
00542     Algo = 5;
00543     break;
00544   }
00545   if ( Algo >= 4 ) {
00546     isFithStep = true;
00547   }
00548 
00549   return isFithStep;
00550 }
00551 // -- method to apply gsf electron selection to EcalDriven seeds
00552 bool 
00553 PFElecTkProducer::applySelection(reco::GsfTrack gsftk) {
00554   if (&(*gsftk.seedRef())==0) return false;
00555   ElectronSeedRef ElSeedRef=gsftk.extra()->seedRef().castTo<ElectronSeedRef>();
00556 
00557   bool passCut = false;
00558   if (ElSeedRef->ctfTrack().isNull()){
00559     if(ElSeedRef->caloCluster().isNull()) return passCut;
00560     SuperClusterRef scRef = ElSeedRef->caloCluster().castTo<SuperClusterRef>();
00561     //do this just to know if exist a SC? 
00562     if(scRef.isNonnull()) {
00563       float caloEne = scRef->energy();
00564       float feta = fabs(scRef->eta()-gsftk.etaMode());
00565       float fphi = fabs(scRef->phi()-gsftk.phiMode());
00566       if (fphi>TMath::Pi()) fphi-= TMath::TwoPi();
00567       if(caloEne > SCEne_ && feta < detaGsfSC_ && fabs(fphi) < dphiGsfSC_)
00568         passCut = true;
00569     }
00570   }
00571   else {
00572     // get all the gsf found by tracker driven
00573     passCut = true;
00574   }
00575   return passCut;
00576 }
00577 bool 
00578 PFElecTkProducer::resolveGsfTracks(const vector<reco::GsfPFRecTrack>  & GsfPFVec, 
00579                                    unsigned int ngsf, 
00580                                    vector<unsigned int> &secondaries,
00581                                    const reco::PFClusterCollection & theEClus) {
00582   bool debugCleaning = debugGsfCleaning_;
00583   bool n_keepGsf = true;
00584 
00585   reco::GsfTrackRef nGsfTrack = GsfPFVec[ngsf].gsfTrackRef();
00586   
00587   if (&(*nGsfTrack->seedRef())==0) return false;    
00588   ElectronSeedRef nElSeedRef=nGsfTrack->extra()->seedRef().castTo<ElectronSeedRef>();
00589 
00590 
00591   TrajectoryStateOnSurface inTSOS = mtsTransform_.innerStateOnSurface((*nGsfTrack));
00592   GlobalVector ninnMom;
00593   float nPin =  nGsfTrack->pMode();
00594   if(inTSOS.isValid()){
00595     mtsMode_->momentumFromModeCartesian(inTSOS,ninnMom);
00596     nPin = ninnMom.mag();
00597   }
00598 
00599   float neta = nGsfTrack->innerMomentum().eta();
00600   float nphi = nGsfTrack->innerMomentum().phi();
00601   
00602 
00603 
00604   
00605   if(debugCleaning)
00606     cout << " PFElecTkProducer:: considering track " << nGsfTrack->pt() 
00607          << " eta,phi " <<  nGsfTrack->eta() << ", " <<  nGsfTrack->phi()  << endl;
00608   
00609   
00610   for (unsigned int igsf=0; igsf< GsfPFVec.size();igsf++) {
00611     if(igsf != ngsf ) {
00612       reco::GsfTrackRef iGsfTrack = GsfPFVec[igsf].gsfTrackRef();
00613 
00614       if(debugCleaning)
00615         cout << " PFElecTkProducer:: and  comparing with track " << iGsfTrack->pt() 
00616              << " eta,phi " <<  iGsfTrack->eta() << ", " <<  iGsfTrack->phi()  << endl;
00617       
00618       float ieta = iGsfTrack->innerMomentum().eta();
00619       float iphi = iGsfTrack->innerMomentum().phi();
00620       float feta = fabs(neta - ieta);
00621       float fphi = fabs(nphi - iphi);
00622       if (fphi>TMath::Pi()) fphi-= TMath::TwoPi();     
00623   
00624 
00625       // apply a superloose preselection only to avoid un-useful cpu time: hard-coded for this reason
00626       if(feta < 0.5 && fabs(fphi) < 1.0) {
00627         if(debugCleaning)
00628           cout << " Entering angular superloose preselection " << endl;
00629         
00630         TrajectoryStateOnSurface i_inTSOS = mtsTransform_.innerStateOnSurface((*iGsfTrack));
00631         GlobalVector i_innMom;
00632         float iPin = iGsfTrack->pMode();
00633         if(i_inTSOS.isValid()){
00634           mtsMode_->momentumFromModeCartesian(i_inTSOS,i_innMom);  
00635           iPin = i_innMom.mag();
00636         }
00637 
00638         if (&(*iGsfTrack->seedRef())==0) continue;   
00639         ElectronSeedRef iElSeedRef=iGsfTrack->extra()->seedRef().castTo<ElectronSeedRef>();
00640 
00641         float SCEnergy = -1.;
00642         // Check if two tracks match the same SC     
00643         bool areBothGsfEcalDriven = false;;
00644         bool isSameSC = isSameEgSC(nElSeedRef,iElSeedRef,areBothGsfEcalDriven,SCEnergy);
00645         
00646         // CASE1 both GsfTracks ecalDriven and match the same SC
00647         if(areBothGsfEcalDriven ) {
00648           if(isSameSC) {
00649             float nEP = SCEnergy/nPin;
00650             float iEP =  SCEnergy/iPin;
00651             if(debugCleaning)
00652               cout << " Entering SAME supercluster case " 
00653                    << " nEP " << nEP 
00654                    << " iEP " << iEP << endl;
00655             
00656             
00657             
00658             // if same SC take the closest or if same
00659             // distance the best E/p
00660             
00661             // Innermost using LostHits technology
00662             bool isSameLayer = false;
00663             bool iGsfInnermostWithLostHits = 
00664               isInnerMostWithLostHits(nGsfTrack,iGsfTrack,isSameLayer);
00665             
00666             
00667             if(debugCleaning)
00668               cout << " iGsf is InnerMostWithLostHits " << iGsfInnermostWithLostHits
00669                    << " isSameLayer " << isSameLayer  << endl;
00670             
00671             if (iGsfInnermostWithLostHits) {
00672               n_keepGsf = false;
00673               return n_keepGsf;
00674             }
00675             else if(isSameLayer){
00676               if(fabs(iEP-1) < fabs(nEP-1)) {
00677                 n_keepGsf = false;
00678                 return n_keepGsf;
00679               }
00680               else{
00681                 secondaries.push_back(igsf);
00682               }
00683             }
00684             else {
00685               // save secondaries gsf track (put selection)
00686             secondaries.push_back(igsf);
00687             }
00688           } // end same SC case
00689         }
00690         else {
00691           // enter in the condition where at least one track is trackerDriven
00692           float minBremDphi =  minTangDist(GsfPFVec[ngsf],GsfPFVec[igsf]);
00693           float nETot = 0.;
00694           float iETot = 0.;
00695           bool isBothGsfTrackerDriven = false;
00696           bool nEcalDriven = false;
00697           bool iEcalDriven = false;
00698           bool isSameScEgPf = isSharingEcalEnergyWithEgSC(GsfPFVec[ngsf],
00699                                                           GsfPFVec[igsf],
00700                                                           nElSeedRef,
00701                                                           iElSeedRef,
00702                                                           theEClus,
00703                                                           isBothGsfTrackerDriven,
00704                                                           nEcalDriven,
00705                                                           iEcalDriven,
00706                                                           nETot,
00707                                                           iETot);
00708 
00709           // check if the first hit of iGsfTrack < nGsfTrack          
00710           bool isSameLayer = false;
00711           bool iGsfInnermostWithLostHits = 
00712             isInnerMostWithLostHits(nGsfTrack,iGsfTrack,isSameLayer);
00713 
00714           if(isSameScEgPf) {
00715             // CASE 2 : One Gsf has reference to a SC and the other one not or both not
00716                     
00717             if(debugCleaning) {
00718               cout << " Sharing ECAL energy passed " 
00719                    << " nEtot " << nETot 
00720                    << " iEtot " << iETot << endl;
00721               if(isBothGsfTrackerDriven) 
00722                 cout << " Both Track are trackerDriven " << endl;
00723             }
00724 
00725             // Innermost using LostHits technology
00726             if (iGsfInnermostWithLostHits) {
00727               n_keepGsf = false;
00728               return n_keepGsf;
00729             }
00730             else if(isSameLayer){
00731               // Thirt Case:  One Gsf has reference to a SC and the other one not or both not
00732               // gsf tracks starts from the same layer
00733               // check number of sharing modules (at least 50%)
00734               // check number of sharing hits (at least 2)
00735               // check charge flip inner/outer
00736 
00737              
00738                 // they share energy
00739               if(isBothGsfTrackerDriven == false) {
00740                 // if at least one Gsf track is EcalDriven choose that one.
00741                 if(iEcalDriven) {
00742                   n_keepGsf = false;
00743                   return n_keepGsf;
00744                 }
00745                 else {
00746                   secondaries.push_back(igsf);
00747                 }
00748               }
00749               else {
00750                 // if both tracks are tracker driven choose that one with the best E/p
00751                 // with ETot = max(En,Ei)
00752 
00753                 float ETot = -1;
00754                 if(nETot != iETot) {
00755                   if(nETot > iETot)
00756                     ETot = nETot;
00757                   else 
00758                     ETot = iETot;
00759                 }
00760                 else {
00761                   ETot = nETot;
00762                 }
00763                 float nEP = ETot/nPin;
00764                 float iEP = ETot/iPin;
00765                 
00766                 
00767                 if(debugCleaning) 
00768                   cout << " nETot " << nETot
00769                        << " iETot " << iETot 
00770                        << " ETot " << ETot << endl 
00771                        << " nPin " << nPin
00772                        << " iPin " << iPin 
00773                        << " nEP " << nEP 
00774                        << " iEP " << iEP << endl;
00775                 
00776         
00777                 if(fabs(iEP-1) < fabs(nEP-1)) {
00778                   n_keepGsf = false;
00779                   return n_keepGsf;
00780                 }
00781                 else{
00782                   secondaries.push_back(igsf);
00783                 }
00784               }
00785             }
00786             else {
00787               secondaries.push_back(igsf);
00788             }
00789           }
00790           else if(feta < detaCutGsfClean_ && minBremDphi < dphiCutGsfClean_) {
00791             // very close tracks
00792             bool secPushedBack = false;
00793             if(nEcalDriven == false && nETot == 0.) {
00794               n_keepGsf = false;
00795               return n_keepGsf;
00796             }
00797             else if(iEcalDriven == false && iETot == 0.) {
00798               secondaries.push_back(igsf);
00799               secPushedBack = true;
00800             }
00801             if(debugCleaning)
00802               cout << " Close Tracks " 
00803                    << " feta " << feta << " fabs(fphi) " << fabs(fphi) 
00804                    << " minBremDphi " <<  minBremDphi 
00805                    << " nETot " << nETot 
00806                    << " iETot " << iETot 
00807                    << " nLostHits " <<  nGsfTrack->trackerExpectedHitsInner().numberOfLostHits() 
00808                    << " iLostHits " << iGsfTrack->trackerExpectedHitsInner().numberOfLostHits() << endl;
00809             
00810             // apply selection only if one track has lost hits
00811             if(applyAngularGsfClean_) {
00812               if (iGsfInnermostWithLostHits) {
00813                 n_keepGsf = false;
00814                 return n_keepGsf;
00815               }
00816               else if(isSameLayer == false) {
00817                 if(secPushedBack == false) 
00818                   secondaries.push_back(igsf);
00819               }
00820             }
00821           }
00822           else if(feta < 0.1 && minBremDphi < 0.2){
00823             // failed all the conditions, discard only tracker driven tracks
00824             // with no PFClusters linked. 
00825             if(debugCleaning)
00826               cout << " Close Tracks and failed all the conditions " 
00827                    << " feta " << feta << " fabs(fphi) " << fabs(fphi) 
00828                    << " minBremDphi " <<  minBremDphi 
00829                    << " nETot " << nETot 
00830                    << " iETot " << iETot 
00831                    << " nLostHits " <<  nGsfTrack->trackerExpectedHitsInner().numberOfLostHits() 
00832                    << " iLostHits " << iGsfTrack->trackerExpectedHitsInner().numberOfLostHits() << endl;
00833             
00834             if(nEcalDriven == false && nETot == 0.) {
00835               n_keepGsf = false;
00836               return n_keepGsf;
00837             }
00838             // Here I do not push back the secondary because considered fakes...
00839           }
00840         }
00841       }
00842     }
00843   }
00844   
00845   return n_keepGsf;
00846 }
00847 float 
00848 PFElecTkProducer::minTangDist(const reco::GsfPFRecTrack& primGsf,
00849                               const reco::GsfPFRecTrack& secGsf) {
00850 
00851   float minDphi = 1000.; 
00852 
00853 
00854   std::vector<reco::PFBrem> primPFBrem = primGsf.PFRecBrem();
00855   std::vector<reco::PFBrem> secPFBrem = secGsf.PFRecBrem();
00856 
00857 
00858   unsigned int cbrem = 0;
00859   for (unsigned isbrem = 0; isbrem < secPFBrem.size(); isbrem++) {
00860     if(secPFBrem[isbrem].indTrajPoint() == 99) continue;
00861     const reco::PFTrajectoryPoint& atSecECAL 
00862       = secPFBrem[isbrem].extrapolatedPoint( reco::PFTrajectoryPoint::ECALEntrance );
00863     if( ! atSecECAL.isValid() ) continue;
00864     float secPhi  = atSecECAL.positionREP().Phi();
00865 
00866     unsigned int sbrem = 0;
00867     for(unsigned ipbrem = 0; ipbrem < primPFBrem.size(); ipbrem++) {
00868       if(primPFBrem[ipbrem].indTrajPoint() == 99) continue;
00869       const reco::PFTrajectoryPoint& atPrimECAL 
00870         = primPFBrem[ipbrem].extrapolatedPoint( reco::PFTrajectoryPoint::ECALEntrance );
00871       if( ! atPrimECAL.isValid() ) continue;
00872       sbrem++;
00873       if(sbrem <= 3) {
00874         float primPhi = atPrimECAL.positionREP().Phi();
00875         
00876         float dphi = fabs(primPhi - secPhi);
00877         if (dphi>TMath::Pi()) dphi-= TMath::TwoPi();     
00878         if(fabs(dphi) < minDphi) {         
00879           minDphi = fabs(dphi);
00880         }
00881       }
00882     }
00883     
00884     
00885     cbrem++;
00886     if(cbrem == 3) 
00887       break;
00888   }
00889   return minDphi;
00890 }
00891 bool 
00892 PFElecTkProducer::isSameEgSC(const reco::ElectronSeedRef& nSeedRef,
00893                              const reco::ElectronSeedRef& iSeedRef,
00894                              bool& bothGsfEcalDriven,
00895                              float& SCEnergy) {
00896   
00897   bool isSameSC = false;
00898 
00899   if(nSeedRef->caloCluster().isNonnull() && iSeedRef->caloCluster().isNonnull()) {
00900     SuperClusterRef nscRef = nSeedRef->caloCluster().castTo<SuperClusterRef>();
00901     SuperClusterRef iscRef = iSeedRef->caloCluster().castTo<SuperClusterRef>();
00902 
00903     if(nscRef.isNonnull() && iscRef.isNonnull()) {
00904       bothGsfEcalDriven = true;
00905       if(nscRef == iscRef) {
00906         isSameSC = true;
00907         // retrieve the supercluster energy
00908         SCEnergy = nscRef->energy();
00909       }
00910     }
00911   }
00912   return isSameSC;
00913 }
00914 bool 
00915 PFElecTkProducer::isSharingEcalEnergyWithEgSC(const reco::GsfPFRecTrack& nGsfPFRecTrack,
00916                                               const reco::GsfPFRecTrack& iGsfPFRecTrack,
00917                                               const reco::ElectronSeedRef& nSeedRef,
00918                                               const reco::ElectronSeedRef& iSeedRef,
00919                                               const reco::PFClusterCollection& theEClus,
00920                                               bool& bothGsfTrackerDriven,
00921                                               bool& nEcalDriven,
00922                                               bool& iEcalDriven,
00923                                               float& nEnergy,
00924                                               float& iEnergy) {
00925   
00926   bool isSharingEnergy = false;
00927 
00928   //which is EcalDriven?
00929   bool oneEcalDriven = true;
00930   SuperClusterRef scRef;
00931   GsfPFRecTrack gsfPfTrack;
00932 
00933   if(nSeedRef->caloCluster().isNonnull()) {
00934     scRef = nSeedRef->caloCluster().castTo<SuperClusterRef>();
00935     nEnergy = scRef->energy();
00936     nEcalDriven = true;
00937     gsfPfTrack = iGsfPFRecTrack;
00938   }
00939   else if(iSeedRef->caloCluster().isNonnull()){
00940     scRef = iSeedRef->caloCluster().castTo<SuperClusterRef>();
00941     iEnergy = scRef->energy();
00942     iEcalDriven = true;
00943     gsfPfTrack = nGsfPFRecTrack;
00944   }
00945   else{
00946      oneEcalDriven = false;
00947   }
00948   
00949   if(oneEcalDriven) {
00950     //run a basic reconstruction for the particle flow
00951 
00952     vector<PFCluster> vecPFClusters;
00953     vecPFClusters.clear();
00954 
00955     for (PFClusterCollection::const_iterator clus = theEClus.begin();
00956          clus != theEClus.end();
00957          clus++ ) {
00958       PFCluster clust = *clus;
00959       clust.calculatePositionREP();
00960 
00961       float deta = fabs(scRef->position().eta() - clust.position().eta());
00962       float dphi = fabs(scRef->position().phi() - clust.position().phi());
00963       if (dphi>TMath::Pi()) dphi-= TMath::TwoPi();
00964  
00965       // Angle preselection between the supercluster and pfclusters
00966       // this is needed just to save some cpu-time for this is hard-coded     
00967       if(deta < 0.5 && fabs(dphi) < 1.0) {
00968         bool foundLink = false;
00969         double distGsf = gsfPfTrack.extrapolatedPoint( reco::PFTrajectoryPoint::ECALShowerMax ).isValid() ?
00970           LinkByRecHit::testTrackAndClusterByRecHit(gsfPfTrack , clust ) : -1.;
00971         // check if it touch the GsfTrack
00972         if(distGsf > 0.) {
00973           if(nEcalDriven) 
00974             iEnergy += clust.energy();
00975           else
00976             nEnergy += clust.energy();
00977           vecPFClusters.push_back(clust);
00978           foundLink = true;
00979         }
00980         // check if it touch the Brem-tangents
00981         if(foundLink == false) {
00982           vector<PFBrem> primPFBrem = gsfPfTrack.PFRecBrem();
00983           for(unsigned ipbrem = 0; ipbrem < primPFBrem.size(); ipbrem++) {
00984             if(primPFBrem[ipbrem].indTrajPoint() == 99) continue;
00985             const reco::PFRecTrack& pfBremTrack = primPFBrem[ipbrem];
00986             double dist = pfBremTrack.extrapolatedPoint( reco::PFTrajectoryPoint::ECALShowerMax ).isValid() ?
00987               LinkByRecHit::testTrackAndClusterByRecHit(pfBremTrack , clust, true ) : -1.;
00988             if(dist > 0.) {
00989               if(nEcalDriven) 
00990                 iEnergy += clust.energy();
00991               else
00992                 nEnergy += clust.energy();
00993               vecPFClusters.push_back(clust);
00994               foundLink = true;
00995             }
00996           }     
00997         }  
00998       } // END if anble preselection
00999     } // PFClusters Loop
01000     if(vecPFClusters.size() > 0 ) {
01001       for(unsigned int pf = 0; pf < vecPFClusters.size(); pf++) {
01002         bool isCommon = ClusterClusterMapping::overlap(vecPFClusters[pf],*scRef);
01003         if(isCommon) {
01004           isSharingEnergy = true;
01005         }
01006         break;
01007       }
01008     }
01009   }
01010   else {
01011     // both tracks are trackerDriven, try ECAL energy matching also in this case.
01012 
01013     bothGsfTrackerDriven = true;
01014     vector<PFCluster> nPFCluster;
01015     vector<PFCluster> iPFCluster;
01016 
01017     nPFCluster.clear();
01018     iPFCluster.clear();
01019 
01020     for (PFClusterCollection::const_iterator clus = theEClus.begin();
01021          clus != theEClus.end();
01022          clus++ ) {
01023       PFCluster clust = *clus;
01024       clust.calculatePositionREP();
01025       
01026       float ndeta = fabs(nGsfPFRecTrack.gsfTrackRef()->eta() - clust.position().eta());
01027       float ndphi = fabs(nGsfPFRecTrack.gsfTrackRef()->phi() - clust.position().phi());
01028       if (ndphi>TMath::Pi()) ndphi-= TMath::TwoPi();
01029       // Apply loose preselection with the track
01030       // just to save cpu time, for this hard-coded
01031       if(ndeta < 0.5 && fabs(ndphi) < 1.0) {
01032         bool foundNLink = false;
01033         
01034         double distGsf = nGsfPFRecTrack.extrapolatedPoint( reco::PFTrajectoryPoint::ECALShowerMax ).isValid() ?
01035           LinkByRecHit::testTrackAndClusterByRecHit(nGsfPFRecTrack , clust ) : -1.;
01036         if(distGsf > 0.) {
01037           nPFCluster.push_back(clust);
01038           nEnergy += clust.energy();
01039           foundNLink = true;
01040         }
01041         if(foundNLink == false) {
01042           const vector<PFBrem>& primPFBrem = nGsfPFRecTrack.PFRecBrem();
01043           for(unsigned ipbrem = 0; ipbrem < primPFBrem.size(); ipbrem++) {
01044             if(primPFBrem[ipbrem].indTrajPoint() == 99) continue;
01045             const reco::PFRecTrack& pfBremTrack = primPFBrem[ipbrem];
01046             if(foundNLink == false) {
01047               double dist = pfBremTrack.extrapolatedPoint( reco::PFTrajectoryPoint::ECALShowerMax ).isValid() ?
01048                 LinkByRecHit::testTrackAndClusterByRecHit(pfBremTrack , clust, true ) : -1.;
01049               if(dist > 0.) {
01050                 nPFCluster.push_back(clust);
01051                 nEnergy += clust.energy();
01052                 foundNLink = true;
01053               }
01054             }
01055           }
01056         }
01057       }
01058 
01059       float ideta = fabs(iGsfPFRecTrack.gsfTrackRef()->eta() - clust.position().eta());
01060       float idphi = fabs(iGsfPFRecTrack.gsfTrackRef()->phi() - clust.position().phi());
01061       if (idphi>TMath::Pi()) idphi-= TMath::TwoPi();
01062       // Apply loose preselection with the track
01063       // just to save cpu time, for this hard-coded
01064       if(ideta < 0.5 && fabs(idphi) < 1.0) {
01065         bool foundILink = false;
01066         double dist = iGsfPFRecTrack.extrapolatedPoint( reco::PFTrajectoryPoint::ECALShowerMax ).isValid() ?
01067           LinkByRecHit::testTrackAndClusterByRecHit(iGsfPFRecTrack , clust ) : -1.;
01068         if(dist > 0.) {
01069           iPFCluster.push_back(clust);
01070           iEnergy += clust.energy();
01071           foundILink = true;
01072         }
01073         if(foundILink == false) {
01074           vector<PFBrem> primPFBrem = iGsfPFRecTrack.PFRecBrem();
01075           for(unsigned ipbrem = 0; ipbrem < primPFBrem.size(); ipbrem++) {
01076             if(primPFBrem[ipbrem].indTrajPoint() == 99) continue;
01077             const reco::PFRecTrack& pfBremTrack = primPFBrem[ipbrem];
01078             if(foundILink == false) {
01079               double dist = LinkByRecHit::testTrackAndClusterByRecHit(pfBremTrack , clust, true );
01080               if(dist > 0.) {
01081                 iPFCluster.push_back(clust);
01082                 iEnergy += clust.energy();
01083                 foundILink = true;
01084               }
01085             }
01086           }
01087         }
01088       }
01089     }
01090 
01091 
01092     if(nPFCluster.size() > 0 && iPFCluster.size() > 0) {
01093       for(unsigned int npf = 0; npf < nPFCluster.size(); npf++) {
01094         for(unsigned int ipf = 0; ipf < iPFCluster.size(); ipf++) {
01095           bool isCommon = ClusterClusterMapping::overlap(nPFCluster[npf],iPFCluster[ipf]);
01096           if(isCommon) {
01097             isSharingEnergy = true;
01098             break;
01099           }
01100         }
01101         if(isSharingEnergy)
01102           break;
01103       }
01104     }
01105   }
01106 
01107   return isSharingEnergy;
01108 }
01109 bool PFElecTkProducer::isInnerMost(const reco::GsfTrackRef& nGsfTrack,
01110                                    const reco::GsfTrackRef& iGsfTrack,
01111                                    bool& sameLayer) {
01112   
01113   // copied by the class RecoEgamma/EgammaElectronAlgos/src/EgAmbiguityTools.cc
01114   // obsolete but the code is kept: now using lost hits method
01115 
01116   reco::HitPattern gsfHitPattern1 = nGsfTrack->hitPattern();
01117   reco::HitPattern gsfHitPattern2 = iGsfTrack->hitPattern();
01118   
01119   // retrieve first valid hit
01120   int gsfHitCounter1 = 0 ;
01121   trackingRecHit_iterator elHitsIt1 ;
01122   for
01123     ( elHitsIt1 = nGsfTrack->recHitsBegin() ;
01124      elHitsIt1 != nGsfTrack->recHitsEnd() ;
01125      elHitsIt1++, gsfHitCounter1++ )
01126     { if (((**elHitsIt1).isValid())) break ; }
01127   
01128   int gsfHitCounter2 = 0 ;
01129   trackingRecHit_iterator elHitsIt2 ;
01130   for
01131     ( elHitsIt2 = iGsfTrack->recHitsBegin() ;
01132      elHitsIt2 != iGsfTrack->recHitsEnd() ;
01133      elHitsIt2++, gsfHitCounter2++ )
01134     { if (((**elHitsIt2).isValid())) break ; }
01135   
01136   uint32_t gsfHit1 = gsfHitPattern1.getHitPattern(gsfHitCounter1) ;
01137   uint32_t gsfHit2 = gsfHitPattern2.getHitPattern(gsfHitCounter2) ;
01138   
01139   
01140   if (gsfHitPattern1.getSubStructure(gsfHit1)!=gsfHitPattern2.getSubStructure(gsfHit2))
01141    { 
01142      return (gsfHitPattern2.getSubStructure(gsfHit2)<gsfHitPattern1.getSubStructure(gsfHit1)); 
01143    }
01144   else if (gsfHitPattern1.getLayer(gsfHit1)!=gsfHitPattern2.getLayer(gsfHit2))
01145     { 
01146       return (gsfHitPattern2.getLayer(gsfHit2)<gsfHitPattern1.getLayer(gsfHit1)); 
01147     }
01148   else
01149    { 
01150      sameLayer = true;
01151      return  false; 
01152    }
01153 }
01154 bool PFElecTkProducer::isInnerMostWithLostHits(const reco::GsfTrackRef& nGsfTrack,
01155                                                const reco::GsfTrackRef& iGsfTrack,
01156                                                bool& sameLayer) {
01157   
01158   // define closest using the lost hits on the expectedhitsineer
01159   unsigned int nLostHits = nGsfTrack->trackerExpectedHitsInner().numberOfLostHits();
01160   unsigned int iLostHits = iGsfTrack->trackerExpectedHitsInner().numberOfLostHits();
01161   
01162   if (nLostHits!=iLostHits) {
01163     return (nLostHits > iLostHits);
01164   } 
01165   else {
01166     sameLayer = true;
01167     return  false; 
01168   }
01169 }
01170 
01171 
01172 
01173 // ------------ method called once each job just before starting event loop  ------------
01174 void 
01175 PFElecTkProducer::beginRun(edm::Run& run,
01176                            const EventSetup& iSetup)
01177 {
01178   ESHandle<MagneticField> magneticField;
01179   iSetup.get<IdealMagneticFieldRecord>().get(magneticField);
01180 
01181   ESHandle<TrackerGeometry> tracker;
01182   iSetup.get<TrackerDigiGeometryRecord>().get(tracker);
01183 
01184 
01185   mtsTransform_ = MultiTrajectoryStateTransform(tracker.product(),magneticField.product());
01186   
01187 
01188   pfTransformer_= new PFTrackTransformer(math::XYZVector(magneticField->inTesla(GlobalPoint(0,0,0))));
01189   
01190   
01191   edm::ESHandle<TransientTrackBuilder> builder;
01192   iSetup.get<TransientTrackRecord>().get("TransientTrackBuilder", builder);
01193   TransientTrackBuilder thebuilder = *(builder.product());
01194   
01195 
01196   if(useConvBremFinder_) {
01197     FILE * fileConvBremID = fopen(path_mvaWeightFileConvBrem_.c_str(), "r");
01198     if (fileConvBremID) {
01199       fclose(fileConvBremID);
01200     }
01201     else {
01202       string err = "PFElecTkProducer: cannot open weight file '";
01203       err += path_mvaWeightFileConvBrem_;
01204       err += "'";
01205       throw invalid_argument( err );
01206     }
01207   }
01208   convBremFinder_ = new ConvBremPFTrackFinder(thebuilder,mvaConvBremFinderID_,path_mvaWeightFileConvBrem_);
01209 
01210 }
01211 
01212 // ------------ method called once each job just after ending the event loop  ------------
01213 void 
01214 PFElecTkProducer::endRun() {
01215   delete pfTransformer_;
01216 }
01217 
01218 //define this as a plug-in
/data/refman/pasoursint/CMSSW_6_1_2_SLHC2/src/RecoParticleFlow/PFTracking/plugins/PFElecTkProducer.cc
