Reference Manual

Go to the documentation of this file.
00001 
00002 #include "RecoEgamma/EgammaElectronAlgos/interface/GsfElectronAlgo.h"
00003 #include "RecoEgamma/EgammaElectronAlgos/interface/EgAmbiguityTools.h"
00004 #include "RecoEgamma/EgammaElectronAlgos/interface/ElectronClassification.h"
00005 #include "RecoEgamma/EgammaElectronAlgos/interface/ElectronMomentumCorrector.h"
00006 #include "RecoEgamma/EgammaElectronAlgos/interface/ElectronEnergyCorrector.h"
00007 #include "RecoEgamma/EgammaElectronAlgos/interface/ElectronUtilities.h"
00008 #include "RecoEgamma/EgammaTools/interface/ConversionFinder.h"
00009 
00010 #include "RecoEcal/EgammaCoreTools/interface/EcalClusterFunctionFactory.h"
00011 #include "RecoEcal/EgammaCoreTools/interface/EcalClusterTools.h"
00012 
00013 #include "DataFormats/ParticleFlowReco/interface/GsfPFRecTrackFwd.h"
00014 #include "DataFormats/ParticleFlowReco/interface/GsfPFRecTrack.h"
00015 #include "DataFormats/GsfTrackReco/interface/GsfTrackFwd.h"
00016 #include "DataFormats/EgammaReco/interface/BasicCluster.h"
00017 #include "DataFormats/EgammaReco/interface/ElectronSeed.h"
00018 #include "DataFormats/TrackCandidate/interface/TrackCandidate.h"
00019 #include "DataFormats/TrackCandidate/interface/TrackCandidateCollection.h"
00020 #include "DataFormats/Math/interface/LorentzVector.h"
00021 #include "DataFormats/TrackingRecHit/interface/TrackingRecHitFwd.h"
00022 #include "DataFormats/EcalDetId/interface/EcalSubdetector.h"
00023 #include "DataFormats/CaloTowers/interface/CaloTowerCollection.h"
00024 #include "DataFormats/EcalDetId/interface/EBDetId.h"
00025 #include "DataFormats/EcalDetId/interface/EEDetId.h"
00026 #include "DataFormats/EcalRecHit/interface/EcalRecHit.h"
00027 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
00028 #include "DataFormats/GeometryVector/interface/GlobalVector.h"
00029 #include "DataFormats/Provenance/interface/Provenance.h"
00030 
00031 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
00032 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
00033 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
00034 #include "Geometry/Records/interface/CaloGeometryRecord.h"
00035 #include "Geometry/Records/interface/CaloTopologyRecord.h"
00036 
00037 #include "TrackingTools/PatternTools/interface/Trajectory.h"
00038 #include "TrackingTools/Records/interface/TrackingComponentsRecord.h"
00039 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
00040 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
00041 #include "TrackingTools/GsfTools/interface/MultiTrajectoryStateTransform.h"
00042 #include "TrackingTools/GsfTools/interface/MultiTrajectoryStateMode.h"
00043 #include "TrackingTools/TransientTrackingRecHit/interface/RecHitComparatorByPosition.h"
00044 
00045 #include "RecoTracker/Record/interface/TrackerRecoGeometryRecord.h"
00046 #include "RecoTracker/Record/interface/CkfComponentsRecord.h"
00047 
00048 #include "CondFormats/DataRecord/interface/EcalChannelStatusRcd.h"
00049 
00050 #include "Geometry/CommonDetUnit/interface/TrackingGeometry.h"
00051 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
00052 
00053 
00054 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00055 
00056 
00057 #include "CLHEP/Units/GlobalPhysicalConstants.h"
00058 
00059 #include <TMath.h>
00060 #include <Math/VectorUtil.h>
00061 #include <Math/Point3D.h>
00062 #include <sstream>
00063 #include <algorithm>
00064 
00065 
00066 using namespace edm ;
00067 using namespace std ;
00068 using namespace reco ;
00069 
00070 
00071 //===================================================================
00072 // GsfElectronAlgo::GeneralData
00073 //===================================================================
00074 
00075 // general data and helpers
00076 struct GsfElectronAlgo::GeneralData
00077  {
00078   // constructors
00079   GeneralData
00080    ( const InputTagsConfiguration &,
00081      const StrategyConfiguration &,
00082      const CutsConfiguration & cutsCfg,
00083      const CutsConfiguration & cutsCfgPflow,
00084      const ElectronHcalHelper::Configuration & hcalCfg,
00085      const ElectronHcalHelper::Configuration & hcalCfgPflow,
00086      const IsolationConfiguration &,
00087      const SpikeConfiguration &,
00088      EcalClusterFunctionBaseClass * superClusterErrorFunction ) ;
00089   ~GeneralData() ;
00090 
00091   // configurables
00092   const InputTagsConfiguration inputCfg ;
00093   const StrategyConfiguration strategyCfg ;
00094   const CutsConfiguration cutsCfg ;
00095   const CutsConfiguration cutsCfgPflow ;
00096   const IsolationConfiguration isoCfg ;
00097   const SpikeConfiguration spikeCfg ;
00098 
00099   // additional configuration and helpers
00100   ElectronHcalHelper * hcalHelper, * hcalHelperPflow ;
00101   EcalClusterFunctionBaseClass * superClusterErrorFunction ;
00102  } ;
00103 
00104  GsfElectronAlgo::GeneralData::GeneralData
00105  ( const InputTagsConfiguration & inputConfig,
00106    const StrategyConfiguration & strategyConfig,
00107    const CutsConfiguration & cutsConfig,
00108    const CutsConfiguration & cutsConfigPflow,
00109    const ElectronHcalHelper::Configuration & hcalConfig,
00110    const ElectronHcalHelper::Configuration & hcalConfigPflow,
00111    const IsolationConfiguration & isoConfig,
00112    const SpikeConfiguration & spikeConfig,
00113    EcalClusterFunctionBaseClass * superClusterErrorFunc
00114  )
00115  : inputCfg(inputConfig),
00116    strategyCfg(strategyConfig),
00117    cutsCfg(cutsConfig),
00118    cutsCfgPflow(cutsConfigPflow),
00119    isoCfg(isoConfig),
00120    spikeCfg(spikeConfig),
00121    hcalHelper(new ElectronHcalHelper(hcalConfig)),
00122    hcalHelperPflow(new ElectronHcalHelper(hcalConfigPflow)),
00123    superClusterErrorFunction(superClusterErrorFunc)
00124  {}
00125 
00126 GsfElectronAlgo::GeneralData::~GeneralData()
00127  {
00128   delete hcalHelper ;
00129   delete hcalHelperPflow ;
00130  }
00131 
00132 //===================================================================
00133 // GsfElectronAlgo::EventSetupData
00134 //===================================================================
00135 
00136 struct GsfElectronAlgo::EventSetupData
00137  {
00138    EventSetupData() ;
00139    ~EventSetupData() ;
00140 
00141    unsigned long long cacheIDGeom ;
00142    unsigned long long cacheIDTopo ;
00143    unsigned long long cacheIDTDGeom ;
00144    unsigned long long cacheIDMagField ;
00145    unsigned long long cacheChStatus ;
00146    unsigned long long cacheSevLevel ;
00147 
00148    edm::ESHandle<MagneticField> magField ;
00149    edm::ESHandle<CaloGeometry> caloGeom ;
00150    edm::ESHandle<CaloTopology> caloTopo ;
00151    edm::ESHandle<TrackerGeometry> trackerHandle ;
00152    edm::ESHandle<EcalChannelStatus> chStatus ;
00153    edm::ESHandle<EcalSeverityLevelAlgo> sevLevel;
00154 
00155    const MultiTrajectoryStateTransform * mtsTransform ;
00156    GsfConstraintAtVertex * constraintAtVtx ;
00157    const MultiTrajectoryStateMode * mtsMode ;
00158 } ;
00159 
00160 GsfElectronAlgo::EventSetupData::EventSetupData()
00161  : cacheIDGeom(0), cacheIDTopo(0), cacheIDTDGeom(0), cacheIDMagField(0),cacheChStatus(0),
00162    cacheSevLevel(0), mtsTransform(0), constraintAtVtx(0), mtsMode(new MultiTrajectoryStateMode)
00163  {}
00164 
00165 GsfElectronAlgo::EventSetupData::~EventSetupData()
00166  {
00167   delete mtsMode ;
00168   delete constraintAtVtx ;
00169   delete mtsTransform ;
00170  }
00171 
00172 
00173 //===================================================================
00174 // GsfElectronAlgo::EventData
00175 //===================================================================
00176 
00177 struct GsfElectronAlgo::EventData
00178  {
00179   // general
00180   edm::Event * event ;
00181   const reco::BeamSpot * beamspot ;
00182   GsfElectronPtrCollection * electrons ;
00183 
00184   EventData() ;
00185   ~EventData() ;
00186 
00187   // utilities
00188   void retreiveOriginalTrackCollections
00189    ( const reco::TrackRef &, const reco::GsfTrackRef & ) ;
00190 
00191   // input collections
00192   edm::Handle<reco::GsfElectronCollection> previousElectrons ;
00193   edm::Handle<reco::GsfElectronCollection> pflowElectrons ;
00194   edm::Handle<reco::GsfElectronCoreCollection> coreElectrons ;
00195   edm::Handle<EcalRecHitCollection> reducedEBRecHits ;
00196   edm::Handle<EcalRecHitCollection> reducedEERecHits ;
00197   edm::Handle<reco::TrackCollection> currentCtfTracks ;
00198   edm::Handle<CaloTowerCollection> towers ;
00199   edm::Handle<edm::ValueMap<float> > pfMva ;
00200   edm::Handle<reco::ElectronSeedCollection> seeds ;
00201   edm::Handle<reco::GsfPFRecTrackCollection> gsfPfRecTracks ;
00202   bool originalCtfTrackCollectionRetreived ;
00203   bool originalGsfTrackCollectionRetreived ;
00204   edm::Handle<reco::TrackCollection> originalCtfTracks ;
00205   edm::Handle<reco::GsfTrackCollection> originalGsfTracks ;
00206 
00207   // isolation helpers
00208   ElectronTkIsolation * tkIsolation03, * tkIsolation04 ;
00209   EgammaTowerIsolation * hadDepth1Isolation03, * hadDepth1Isolation04 ;
00210   EgammaTowerIsolation * hadDepth2Isolation03, * hadDepth2Isolation04 ;
00211   EcalRecHitMetaCollection * ecalBarrelHitsMeta ;
00212   EcalRecHitMetaCollection * ecalEndcapHitsMeta ;
00213   EgammaRecHitIsolation * ecalBarrelIsol03, * ecalBarrelIsol04 ;
00214   EgammaRecHitIsolation * ecalEndcapIsol03, * ecalEndcapIsol04 ;
00215  } ;
00216 
00217 GsfElectronAlgo::EventData::EventData()
00218  : event(0), beamspot(0),
00219    originalCtfTrackCollectionRetreived(false),
00220    originalGsfTrackCollectionRetreived(false),
00221    tkIsolation03(0), tkIsolation04(0),
00222    hadDepth1Isolation03(0), hadDepth1Isolation04(0),
00223    hadDepth2Isolation03(0), hadDepth2Isolation04(0),
00224    ecalBarrelHitsMeta(0), ecalEndcapHitsMeta(0),
00225    ecalBarrelIsol03(0), ecalBarrelIsol04(0),
00226    ecalEndcapIsol03(0), ecalEndcapIsol04(0)
00227  {
00228   electrons = new GsfElectronPtrCollection ;
00229  }
00230 
00231 GsfElectronAlgo::EventData::~EventData()
00232  {
00233   delete tkIsolation03 ;
00234   delete tkIsolation04 ;
00235   delete hadDepth1Isolation03 ;
00236   delete hadDepth1Isolation04 ;
00237   delete hadDepth2Isolation03 ;
00238   delete hadDepth2Isolation04 ;
00239   delete ecalBarrelHitsMeta ;
00240   delete ecalEndcapHitsMeta ;
00241   delete ecalBarrelIsol03 ;
00242   delete ecalBarrelIsol04 ;
00243   delete ecalEndcapIsol03 ;
00244   delete ecalEndcapIsol04 ;
00245 
00246   GsfElectronPtrCollection::const_iterator it ;
00247   for ( it = electrons->begin() ; it != electrons->end() ; it++ )
00248    { delete (*it) ; }
00249   delete electrons ;
00250  }
00251 
00252 void GsfElectronAlgo::EventData::retreiveOriginalTrackCollections
00253  ( const reco::TrackRef & ctfTrack, const reco::GsfTrackRef & gsfTrack )
00254  {
00255   if ((!originalCtfTrackCollectionRetreived)&&(ctfTrack.isNonnull()))
00256    {
00257     event->get(ctfTrack.id(),originalCtfTracks) ;
00258     originalCtfTrackCollectionRetreived = true ;
00259    }
00260   if ((!originalGsfTrackCollectionRetreived)&&(gsfTrack.isNonnull()))
00261    {
00262     event->get(gsfTrack.id(),originalGsfTracks) ;
00263     originalGsfTrackCollectionRetreived = true ;
00264    }
00265  }
00266 
00267 
00268 //===================================================================
00269 // GsfElectronAlgo::ElectronData
00270 //===================================================================
00271 
00272 struct GsfElectronAlgo::ElectronData
00273  {
00274   // Refs to subproducts
00275   const reco::GsfElectronCoreRef coreRef ;
00276   const reco::GsfTrackRef gsfTrackRef ;
00277   const reco::SuperClusterRef superClusterRef ;
00278   reco::TrackRef ctfTrackRef ;
00279   float shFracInnerHits ;
00280   const reco::BeamSpot beamSpot ;
00281 
00282   // constructors
00283   ElectronData
00284    ( const reco::GsfElectronCoreRef & core,
00285      const reco::BeamSpot & bs ) ;
00286   ~ElectronData() ;
00287 
00288   // utilities
00289   void checkCtfTrack( edm::Handle<reco::TrackCollection> currentCtfTracks ) ;
00290   void computeCharge( int & charge, reco::GsfElectron::ChargeInfo & info ) ;
00291   CaloClusterPtr getEleBasicCluster( const MultiTrajectoryStateTransform * ) ;
00292   bool calculateTSOS( const MultiTrajectoryStateTransform *, GsfConstraintAtVertex * ) ;
00293   void calculateMode( const MultiTrajectoryStateMode * mtsMode ) ;
00294   Candidate::LorentzVector calculateMomentum() ;
00295 
00296   // TSOS
00297   TrajectoryStateOnSurface innTSOS ;
00298   TrajectoryStateOnSurface outTSOS ;
00299   TrajectoryStateOnSurface vtxTSOS ;
00300   TrajectoryStateOnSurface sclTSOS ;
00301   TrajectoryStateOnSurface seedTSOS ;
00302   TrajectoryStateOnSurface eleTSOS ;
00303   TrajectoryStateOnSurface constrainedVtxTSOS ;
00304 
00305   // mode
00306   GlobalVector innMom, seedMom, eleMom, sclMom, vtxMom, outMom ;
00307   GlobalPoint innPos, seedPos, elePos, sclPos, vtxPos, outPos ;
00308   GlobalVector vtxMomWithConstraint ;
00309  } ;
00310 
00311 GsfElectronAlgo::ElectronData::ElectronData
00312  ( const reco::GsfElectronCoreRef & core,
00313    const reco::BeamSpot & bs )
00314  : coreRef(core),
00315    gsfTrackRef(coreRef->gsfTrack()),
00316    superClusterRef(coreRef->superCluster()),
00317    ctfTrackRef(coreRef->ctfTrack()), shFracInnerHits(coreRef->ctfGsfOverlap()),
00318    beamSpot(bs)
00319  {}
00320 
00321 GsfElectronAlgo::ElectronData::~ElectronData()
00322  {}
00323 
00324 void GsfElectronAlgo::ElectronData::checkCtfTrack( edm::Handle<reco::TrackCollection> currentCtfTracks )
00325  {
00326   if (!ctfTrackRef.isNull()) return ;
00327 
00328   // Code below from Puneeth Kalavase
00329 
00330   shFracInnerHits = 0 ;
00331   const TrackCollection * ctfTrackCollection = currentCtfTracks.product() ;
00332 
00333   // get the Hit Pattern for the gsfTrack
00334   const HitPattern & gsfHitPattern = gsfTrackRef->hitPattern() ;
00335 
00336   unsigned int counter ;
00337   TrackCollection::const_iterator ctfTkIter ;
00338   for ( ctfTkIter = ctfTrackCollection->begin() , counter = 0 ;
00339         ctfTkIter != ctfTrackCollection->end() ; ctfTkIter++, counter++ )
00340    {
00341 
00342     double dEta = gsfTrackRef->eta() - ctfTkIter->eta() ;
00343     double dPhi = gsfTrackRef->phi() - ctfTkIter->phi() ;
00344     double pi = acos(-1.);
00345     if(std::abs(dPhi) > pi) dPhi = 2*pi - std::abs(dPhi) ;
00346 
00347     // dont want to look at every single track in the event!
00348     if (sqrt(dEta*dEta + dPhi*dPhi) > 0.3) continue ;
00349 
00350     unsigned int shared = 0 ;
00351     int gsfHitCounter = 0 ;
00352     int numGsfInnerHits = 0 ;
00353     int numCtfInnerHits = 0 ;
00354 
00355     // get the CTF Track Hit Pattern
00356     const HitPattern & ctfHitPattern = ctfTkIter->hitPattern() ;
00357 
00358     trackingRecHit_iterator elHitsIt ;
00359     for ( elHitsIt = gsfTrackRef->recHitsBegin() ;
00360           elHitsIt != gsfTrackRef->recHitsEnd() ;
00361           elHitsIt++, gsfHitCounter++ )
00362      {
00363       if (!((**elHitsIt).isValid()))  //count only valid Hits
00364        { continue ; }
00365 
00366       // look only in the pixels/TIB/TID
00367       uint32_t gsfHit = gsfHitPattern.getHitPattern(gsfHitCounter) ;
00368       if (!(gsfHitPattern.pixelHitFilter(gsfHit) ||
00369           gsfHitPattern.stripTIBHitFilter(gsfHit) ||
00370           gsfHitPattern.stripTIDHitFilter(gsfHit) ) )
00371        { continue ; }
00372 
00373       numGsfInnerHits++ ;
00374 
00375       int ctfHitsCounter = 0 ;
00376       numCtfInnerHits = 0 ;
00377       trackingRecHit_iterator ctfHitsIt ;
00378       for ( ctfHitsIt = ctfTkIter->recHitsBegin() ;
00379             ctfHitsIt != ctfTkIter->recHitsEnd() ;
00380             ctfHitsIt++, ctfHitsCounter++ )
00381        {
00382         if(!((**ctfHitsIt).isValid())) //count only valid Hits!
00383          { continue ; }
00384 
00385         uint32_t ctfHit = ctfHitPattern.getHitPattern(ctfHitsCounter);
00386         if( !(ctfHitPattern.pixelHitFilter(ctfHit) ||
00387               ctfHitPattern.stripTIBHitFilter(ctfHit) ||
00388               ctfHitPattern.stripTIDHitFilter(ctfHit) ) )
00389          { continue ; }
00390 
00391         numCtfInnerHits++ ;
00392 
00393         if( (**elHitsIt).sharesInput(&(**ctfHitsIt),TrackingRecHit::all) )
00394          {
00395           shared++ ;
00396           break ;
00397          }
00398 
00399        } //ctfHits iterator
00400 
00401      } //gsfHits iterator
00402 
00403     if ((numGsfInnerHits==0)||(numCtfInnerHits==0))
00404      { continue ; }
00405 
00406     if ( static_cast<float>(shared)/min(numGsfInnerHits,numCtfInnerHits) > shFracInnerHits )
00407      {
00408       shFracInnerHits = static_cast<float>(shared)/min(numGsfInnerHits, numCtfInnerHits);
00409       ctfTrackRef = TrackRef(currentCtfTracks,counter);
00410      }
00411    } //ctfTrack iterator
00412  }
00413 
00414 void GsfElectronAlgo::ElectronData::computeCharge
00415  ( int & charge, GsfElectron::ChargeInfo & info )
00416  {
00417   // determine charge from SC
00418   GlobalPoint orig, scpos ;
00419   ele_convert(beamSpot.position(),orig) ;
00420   ele_convert(superClusterRef->position(),scpos) ;
00421   GlobalVector scvect(scpos-orig) ;
00422   GlobalPoint inntkpos = innTSOS.globalPosition() ;
00423   GlobalVector inntkvect = GlobalVector(inntkpos-orig) ;
00424   float dPhiInnEle=normalized_phi(scvect.phi()-inntkvect.phi()) ;
00425   if(dPhiInnEle>0) info.scPixCharge = -1 ;
00426   else info.scPixCharge = 1 ;
00427 
00428   // flags
00429   int chargeGsf = gsfTrackRef->charge() ;
00430   info.isGsfScPixConsistent = ((chargeGsf*info.scPixCharge)>0) ;
00431   info.isGsfCtfConsistent = (ctfTrackRef.isNonnull()&&((chargeGsf*ctfTrackRef->charge())>0)) ;
00432   info.isGsfCtfScPixConsistent = (info.isGsfScPixConsistent&&info.isGsfCtfConsistent) ;
00433 
00434   // default charge
00435   if (info.isGsfScPixConsistent||ctfTrackRef.isNull())
00436    { charge = info.scPixCharge ; }
00437   else
00438    { charge = ctfTrackRef->charge() ; }
00439  }
00440 
00441 CaloClusterPtr GsfElectronAlgo::ElectronData::getEleBasicCluster
00442  ( const MultiTrajectoryStateTransform * mtsTransform )
00443  {
00444   CaloClusterPtr eleRef ;
00445   TrajectoryStateOnSurface tempTSOS ;
00446   TrajectoryStateOnSurface outTSOS = mtsTransform->outerStateOnSurface(*gsfTrackRef) ;
00447   float dphimin = 1.e30 ;
00448   for (CaloCluster_iterator bc=superClusterRef->clustersBegin(); bc!=superClusterRef->clustersEnd(); bc++)
00449    {
00450     GlobalPoint posclu((*bc)->position().x(),(*bc)->position().y(),(*bc)->position().z()) ;
00451     tempTSOS = mtsTransform->extrapolatedState(outTSOS,posclu) ;
00452     if (!tempTSOS.isValid()) tempTSOS=outTSOS ;
00453     GlobalPoint extrap = tempTSOS.globalPosition() ;
00454     float dphi = EleRelPointPair(posclu,extrap,beamSpot.position()).dPhi() ;
00455     if (std::abs(dphi)<dphimin)
00456      {
00457       dphimin = std::abs(dphi) ;
00458       eleRef = (*bc);
00459       eleTSOS = tempTSOS ;
00460      }
00461    }
00462   return eleRef ;
00463  }
00464 
00465 bool GsfElectronAlgo::ElectronData::calculateTSOS
00466  ( const MultiTrajectoryStateTransform * mtsTransform, GsfConstraintAtVertex * constraintAtVtx )
00467  {
00468   //at innermost point
00469   innTSOS = mtsTransform->innerStateOnSurface(*gsfTrackRef);
00470   if (!innTSOS.isValid()) return false;
00471 
00472   //at vertex
00473   // innermost state propagation to the beam spot position
00474   GlobalPoint bsPos ;
00475   ele_convert(beamSpot.position(),bsPos) ;
00476   vtxTSOS = mtsTransform->extrapolatedState(innTSOS,bsPos) ;
00477   if (!vtxTSOS.isValid()) vtxTSOS=innTSOS;
00478 
00479   //at seed
00480   outTSOS = mtsTransform->outerStateOnSurface(*gsfTrackRef);
00481   if (!outTSOS.isValid()) return false;
00482 
00483   //    TrajectoryStateOnSurface seedTSOS
00484   seedTSOS = mtsTransform->extrapolatedState(outTSOS,
00485            GlobalPoint(superClusterRef->seed()->position().x(),
00486                superClusterRef->seed()->position().y(),
00487                  superClusterRef->seed()->position().z()));
00488   if (!seedTSOS.isValid()) seedTSOS=outTSOS;
00489 
00490   // at scl
00491   sclTSOS = mtsTransform->extrapolatedState(innTSOS,GlobalPoint(superClusterRef->x(),superClusterRef->y(),superClusterRef->z()));
00492   if (!sclTSOS.isValid()) sclTSOS=outTSOS;
00493 
00494   // constrained momentum
00495   constrainedVtxTSOS = constraintAtVtx->constrainAtBeamSpot(*gsfTrackRef,beamSpot);
00496 
00497   return true ;
00498  }
00499 
00500 void GsfElectronAlgo::ElectronData::calculateMode( const MultiTrajectoryStateMode * mtsMode )
00501  {
00502   mtsMode->momentumFromModeCartesian(innTSOS,innMom) ;
00503   mtsMode->positionFromModeCartesian(innTSOS,innPos) ;
00504   mtsMode->momentumFromModeCartesian(seedTSOS,seedMom) ;
00505   mtsMode->positionFromModeCartesian(seedTSOS,seedPos) ;
00506   mtsMode->momentumFromModeCartesian(eleTSOS,eleMom) ;
00507   mtsMode->positionFromModeCartesian(eleTSOS,elePos) ;
00508   mtsMode->momentumFromModeCartesian(sclTSOS,sclMom) ;
00509   mtsMode->positionFromModeCartesian(sclTSOS,sclPos) ;
00510   mtsMode->momentumFromModeCartesian(vtxTSOS,vtxMom) ;
00511   mtsMode->positionFromModeCartesian(vtxTSOS,vtxPos) ;
00512   mtsMode->momentumFromModeCartesian(outTSOS,outMom);
00513   mtsMode->positionFromModeCartesian(outTSOS,outPos) ;
00514   mtsMode->momentumFromModeCartesian(constrainedVtxTSOS,vtxMomWithConstraint);
00515  }
00516 
00517 Candidate::LorentzVector GsfElectronAlgo::ElectronData::calculateMomentum()
00518  {
00519   double scale = superClusterRef->energy()/vtxMom.mag() ;
00520   return Candidate::LorentzVector
00521    ( vtxMom.x()*scale,vtxMom.y()*scale,vtxMom.z()*scale,
00522      superClusterRef->energy() ) ;
00523  }
00524 
00525 void GsfElectronAlgo::calculateShowerShape( const reco::SuperClusterRef & theClus, bool pflow, reco::GsfElectron::ShowerShape & showerShape )
00526  {
00527   const reco::CaloCluster & seedCluster = *(theClus->seed()) ;
00528   // temporary, till CaloCluster->seed() is made available
00529   DetId seedXtalId = seedCluster.hitsAndFractions()[0].first ;
00530   int detector = seedXtalId.subdetId() ;
00531 
00532   const CaloTopology * topology = eventSetupData_->caloTopo.product() ;
00533   const CaloGeometry * geometry = eventSetupData_->caloGeom.product() ;
00534   const EcalRecHitCollection * reducedRecHits = 0 ;
00535   if (detector==EcalBarrel)
00536    { reducedRecHits = eventData_->reducedEBRecHits.product() ; }
00537   else
00538    { reducedRecHits = eventData_->reducedEERecHits.product() ; }
00539 
00540   std::vector<float> covariances = EcalClusterTools::covariances(seedCluster,reducedRecHits,topology,geometry) ;
00541   std::vector<float> localCovariances = EcalClusterTools::localCovariances(seedCluster,reducedRecHits,topology) ;
00542   showerShape.sigmaEtaEta = sqrt(covariances[0]) ;
00543   showerShape.sigmaIetaIeta = sqrt(localCovariances[0]) ;
00544   showerShape.e1x5 = EcalClusterTools::e1x5(seedCluster,reducedRecHits,topology)  ;
00545   showerShape.e2x5Max = EcalClusterTools::e2x5Max(seedCluster,reducedRecHits,topology)  ;
00546   showerShape.e5x5 = EcalClusterTools::e5x5(seedCluster,reducedRecHits,topology) ;
00547 
00548   if (pflow)
00549    {
00550     showerShape.hcalDepth1OverEcal = generalData_->hcalHelperPflow->hcalESumDepth1(*theClus)/theClus->energy() ;
00551     showerShape.hcalDepth2OverEcal = generalData_->hcalHelperPflow->hcalESumDepth2(*theClus)/theClus->energy() ;
00552    }
00553   else
00554    {
00555     showerShape.hcalDepth1OverEcal = generalData_->hcalHelper->hcalESumDepth1(*theClus)/theClus->energy() ;
00556     showerShape.hcalDepth2OverEcal = generalData_->hcalHelper->hcalESumDepth2(*theClus)/theClus->energy() ;
00557    }
00558  }
00559 
00560 
00561 //===================================================================
00562 // GsfElectronAlgo
00563 //===================================================================
00564 
00565 GsfElectronAlgo::GsfElectronAlgo
00566  ( const InputTagsConfiguration & inputCfg,
00567    const StrategyConfiguration & strategyCfg,
00568    const CutsConfiguration & cutsCfg,
00569    const CutsConfiguration & cutsCfgPflow,
00570    const ElectronHcalHelper::Configuration & hcalCfg,
00571    const ElectronHcalHelper::Configuration & hcalCfgPflow,
00572    const IsolationConfiguration & isoCfg,
00573    const SpikeConfiguration & spikeCfg,
00574    EcalClusterFunctionBaseClass * superClusterErrorFunction
00575  )
00576  : generalData_(new GeneralData(inputCfg,strategyCfg,cutsCfg,cutsCfgPflow,hcalCfg,hcalCfgPflow,isoCfg,spikeCfg,superClusterErrorFunction)),
00577    eventSetupData_(new EventSetupData),
00578    eventData_(0), electronData_(0)
00579  {}
00580 
00581 GsfElectronAlgo::~GsfElectronAlgo()
00582  {
00583   delete generalData_ ;
00584   delete eventSetupData_ ;
00585   delete eventData_ ;
00586   delete electronData_ ;
00587  }
00588 
00589 void GsfElectronAlgo::checkSetup( const edm::EventSetup & es )
00590  {
00591   // get EventSetupRecords if needed
00592   bool updateField(false);
00593   if (eventSetupData_->cacheIDMagField!=es.get<IdealMagneticFieldRecord>().cacheIdentifier()){
00594     updateField = true;
00595     eventSetupData_->cacheIDMagField=es.get<IdealMagneticFieldRecord>().cacheIdentifier();
00596     es.get<IdealMagneticFieldRecord>().get(eventSetupData_->magField);
00597   }
00598 
00599   bool updateGeometry(false);
00600   if (eventSetupData_->cacheIDTDGeom!=es.get<TrackerDigiGeometryRecord>().cacheIdentifier()){
00601     updateGeometry = true;
00602     eventSetupData_->cacheIDTDGeom=es.get<TrackerDigiGeometryRecord>().cacheIdentifier();
00603     es.get<TrackerDigiGeometryRecord>().get(eventSetupData_->trackerHandle);
00604   }
00605 
00606   if ( updateField || updateGeometry ) {
00607     delete eventSetupData_->mtsTransform ;
00608     eventSetupData_->mtsTransform = new MultiTrajectoryStateTransform(eventSetupData_->trackerHandle.product(),eventSetupData_->magField.product());
00609     delete eventSetupData_->constraintAtVtx ;
00610     eventSetupData_->constraintAtVtx = new GsfConstraintAtVertex(es) ;
00611   }
00612 
00613   if (eventSetupData_->cacheIDGeom!=es.get<CaloGeometryRecord>().cacheIdentifier()){
00614     eventSetupData_->cacheIDGeom=es.get<CaloGeometryRecord>().cacheIdentifier();
00615     es.get<CaloGeometryRecord>().get(eventSetupData_->caloGeom);
00616   }
00617 
00618   if (eventSetupData_->cacheIDTopo!=es.get<CaloTopologyRecord>().cacheIdentifier()){
00619     eventSetupData_->cacheIDTopo=es.get<CaloTopologyRecord>().cacheIdentifier();
00620     es.get<CaloTopologyRecord>().get(eventSetupData_->caloTopo);
00621   }
00622 
00623   generalData_->hcalHelper->checkSetup(es) ;
00624   generalData_->hcalHelperPflow->checkSetup(es) ;
00625 
00626   if (generalData_->superClusterErrorFunction)
00627    { generalData_->superClusterErrorFunction->init(es) ; }
00628 
00629 
00630   if(eventSetupData_->cacheChStatus!=es.get<EcalChannelStatusRcd>().cacheIdentifier()){
00631     eventSetupData_->cacheChStatus=es.get<EcalChannelStatusRcd>().cacheIdentifier();
00632     es.get<EcalChannelStatusRcd>().get(eventSetupData_->chStatus);
00633   }
00634 
00635   if(eventSetupData_->cacheSevLevel != es.get<EcalSeverityLevelAlgoRcd>().cacheIdentifier()){
00636     eventSetupData_->cacheSevLevel = es.get<EcalSeverityLevelAlgoRcd>().cacheIdentifier();
00637     es.get<EcalSeverityLevelAlgoRcd>().get(eventSetupData_->sevLevel);
00638   }
00639  }
00640 
00641 void GsfElectronAlgo::copyElectrons( GsfElectronCollection & outEle )
00642  {
00643   GsfElectronPtrCollection::const_iterator it ;
00644   for
00645    ( it = eventData_->electrons->begin() ;
00646      it != eventData_->electrons->end() ;
00647      it++ )
00648    { outEle.push_back(**it) ; }
00649  }
00650 
00651 void GsfElectronAlgo::beginEvent( edm::Event & event )
00652  {
00653   if (eventData_!=0)
00654    { throw cms::Exception("GsfElectronAlgo|InternalError")<<"unexpected event data" ; }
00655   eventData_ = new EventData ;
00656 
00657   // init the handles linked to the current event
00658   eventData_->event = &event ;
00659   event.getByLabel(generalData_->inputCfg.previousGsfElectrons,eventData_->previousElectrons) ;
00660   event.getByLabel(generalData_->inputCfg.pflowGsfElectronsTag,eventData_->pflowElectrons) ;
00661   event.getByLabel(generalData_->inputCfg.gsfElectronCores,eventData_->coreElectrons) ;
00662   event.getByLabel(generalData_->inputCfg.ctfTracks,eventData_->currentCtfTracks) ;
00663   event.getByLabel(generalData_->inputCfg.reducedBarrelRecHitCollection,eventData_->reducedEBRecHits) ;
00664   event.getByLabel(generalData_->inputCfg.reducedEndcapRecHitCollection,eventData_->reducedEERecHits) ;
00665   event.getByLabel(generalData_->inputCfg.hcalTowersTag,eventData_->towers) ;
00666   event.getByLabel(generalData_->inputCfg.pfMVA,eventData_->pfMva) ;
00667   event.getByLabel(generalData_->inputCfg.seedsTag,eventData_->seeds) ;
00668   if (generalData_->strategyCfg.useGsfPfRecTracks)
00669    { event.getByLabel(generalData_->inputCfg.gsfPfRecTracksTag,eventData_->gsfPfRecTracks) ; }
00670 
00671   // get the beamspot from the Event:
00672   edm::Handle<reco::BeamSpot> recoBeamSpotHandle ;
00673   event.getByLabel(generalData_->inputCfg.beamSpotTag,recoBeamSpotHandle) ;
00674   eventData_->beamspot = recoBeamSpotHandle.product() ;
00675 
00676   // prepare access to hcal data
00677   generalData_->hcalHelper->readEvent(event) ;
00678   generalData_->hcalHelperPflow->readEvent(event) ;
00679 
00680   // Isolation algos
00681   float extRadiusSmall=0.3, extRadiusLarge=0.4 ;
00682   float intRadiusBarrel=generalData_->isoCfg.intRadiusBarrelTk, intRadiusEndcap=generalData_->isoCfg.intRadiusEndcapTk, stripBarrel=generalData_->isoCfg.stripBarrelTk, stripEndcap=generalData_->isoCfg.stripEndcapTk ;
00683   float ptMin=generalData_->isoCfg.ptMinTk, maxVtxDist=generalData_->isoCfg.maxVtxDistTk, drb=generalData_->isoCfg.maxDrbTk;
00684   eventData_->tkIsolation03 = new ElectronTkIsolation(extRadiusSmall,intRadiusBarrel,intRadiusEndcap,stripBarrel,stripEndcap,ptMin,maxVtxDist,drb,eventData_->currentCtfTracks.product(),eventData_->beamspot->position()) ;
00685   eventData_->tkIsolation04 = new ElectronTkIsolation(extRadiusLarge,intRadiusBarrel,intRadiusEndcap,stripBarrel,stripEndcap,ptMin,maxVtxDist,drb,eventData_->currentCtfTracks.product(),eventData_->beamspot->position()) ;
00686 
00687   float egHcalIsoConeSizeOutSmall=0.3, egHcalIsoConeSizeOutLarge=0.4;
00688   float egHcalIsoConeSizeIn=generalData_->isoCfg.intRadiusHcal,egHcalIsoPtMin=generalData_->isoCfg.etMinHcal;
00689   int egHcalDepth1=1, egHcalDepth2=2;
00690   eventData_->hadDepth1Isolation03 = new EgammaTowerIsolation(egHcalIsoConeSizeOutSmall,egHcalIsoConeSizeIn,egHcalIsoPtMin,egHcalDepth1,eventData_->towers.product()) ;
00691   eventData_->hadDepth2Isolation03 = new EgammaTowerIsolation(egHcalIsoConeSizeOutSmall,egHcalIsoConeSizeIn,egHcalIsoPtMin,egHcalDepth2,eventData_->towers.product()) ;
00692   eventData_->hadDepth1Isolation04 = new EgammaTowerIsolation(egHcalIsoConeSizeOutLarge,egHcalIsoConeSizeIn,egHcalIsoPtMin,egHcalDepth1,eventData_->towers.product()) ;
00693   eventData_->hadDepth2Isolation04 = new EgammaTowerIsolation(egHcalIsoConeSizeOutLarge,egHcalIsoConeSizeIn,egHcalIsoPtMin,egHcalDepth2,eventData_->towers.product()) ;
00694 
00695   float egIsoConeSizeOutSmall=0.3, egIsoConeSizeOutLarge=0.4, egIsoJurassicWidth=generalData_->isoCfg.jurassicWidth;
00696   float egIsoPtMinBarrel=generalData_->isoCfg.etMinBarrel,egIsoEMinBarrel=generalData_->isoCfg.eMinBarrel, egIsoConeSizeInBarrel=generalData_->isoCfg.intRadiusEcalBarrel;
00697   float egIsoPtMinEndcap=generalData_->isoCfg.etMinEndcaps,egIsoEMinEndcap=generalData_->isoCfg.eMinEndcaps, egIsoConeSizeInEndcap=generalData_->isoCfg.intRadiusEcalEndcaps;
00698   eventData_->ecalBarrelHitsMeta = new EcalRecHitMetaCollection(*eventData_->reducedEBRecHits) ;
00699   eventData_->ecalEndcapHitsMeta = new EcalRecHitMetaCollection(*eventData_->reducedEERecHits) ;
00700   eventData_->ecalBarrelIsol03 = new EgammaRecHitIsolation(egIsoConeSizeOutSmall,egIsoConeSizeInBarrel,egIsoJurassicWidth,egIsoPtMinBarrel,egIsoEMinBarrel,eventSetupData_->caloGeom,eventData_->ecalBarrelHitsMeta,eventSetupData_->sevLevel.product(),DetId::Ecal);
00701   eventData_->ecalBarrelIsol04 = new EgammaRecHitIsolation(egIsoConeSizeOutLarge,egIsoConeSizeInBarrel,egIsoJurassicWidth,egIsoPtMinBarrel,egIsoEMinBarrel,eventSetupData_->caloGeom,eventData_->ecalBarrelHitsMeta,eventSetupData_->sevLevel.product(),DetId::Ecal);
00702   eventData_->ecalEndcapIsol03 = new EgammaRecHitIsolation(egIsoConeSizeOutSmall,egIsoConeSizeInEndcap,egIsoJurassicWidth,egIsoPtMinEndcap,egIsoEMinEndcap,eventSetupData_->caloGeom,eventData_->ecalEndcapHitsMeta,eventSetupData_->sevLevel.product(),DetId::Ecal);
00703   eventData_->ecalEndcapIsol04 = new EgammaRecHitIsolation(egIsoConeSizeOutLarge,egIsoConeSizeInEndcap,egIsoJurassicWidth,egIsoPtMinEndcap,egIsoEMinEndcap,eventSetupData_->caloGeom,eventData_->ecalEndcapHitsMeta,eventSetupData_->sevLevel.product(),DetId::Ecal);
00704   eventData_->ecalBarrelIsol03->setUseNumCrystals(generalData_->isoCfg.useNumCrystals);
00705   eventData_->ecalBarrelIsol03->setVetoClustered(generalData_->isoCfg.vetoClustered);
00706   //eventData_->ecalBarrelIsol03->doSpikeRemoval(eventData_->reducedEBRecHits.product(),eventSetupData_->chStatus.product(),generalData_->spikeCfg.severityLevelCut,generalData_->spikeCfg.severityRecHitThreshold,generalData_->spikeCfg.spikeId,generalData_->spikeCfg.spikeIdThreshold);
00707   eventData_->ecalBarrelIsol03->doSpikeRemoval(eventData_->reducedEBRecHits.product(),eventSetupData_->chStatus.product(),generalData_->spikeCfg.severityLevelCut);
00708   eventData_->ecalBarrelIsol03->doFlagChecks(generalData_->spikeCfg.recHitFlagsToBeExcluded);
00709   eventData_->ecalBarrelIsol04->setUseNumCrystals(generalData_->isoCfg.useNumCrystals);
00710   eventData_->ecalBarrelIsol04->setVetoClustered(generalData_->isoCfg.vetoClustered);
00711   //eventData_->ecalBarrelIsol04->doSpikeRemoval(eventData_->reducedEBRecHits.product(),eventSetupData_->chStatus.product(),generalData_->spikeCfg.severityLevelCut,generalData_->spikeCfg.severityRecHitThreshold,generalData_->spikeCfg.spikeId,generalData_->spikeCfg.spikeIdThreshold);
00712   eventData_->ecalBarrelIsol04->doSpikeRemoval(eventData_->reducedEBRecHits.product(),eventSetupData_->chStatus.product(),generalData_->spikeCfg.severityLevelCut);
00713   eventData_->ecalBarrelIsol04->doFlagChecks(generalData_->spikeCfg.recHitFlagsToBeExcluded);
00714   eventData_->ecalEndcapIsol03->setUseNumCrystals(generalData_->isoCfg.useNumCrystals);
00715   eventData_->ecalEndcapIsol03->setVetoClustered(generalData_->isoCfg.vetoClustered);
00716   eventData_->ecalEndcapIsol03->doFlagChecks(generalData_->spikeCfg.recHitFlagsToBeExcluded);
00717   eventData_->ecalEndcapIsol04->setUseNumCrystals(generalData_->isoCfg.useNumCrystals);
00718   eventData_->ecalEndcapIsol04->setVetoClustered(generalData_->isoCfg.vetoClustered);
00719   eventData_->ecalEndcapIsol04->doFlagChecks(generalData_->spikeCfg.recHitFlagsToBeExcluded);
00720  }
00721 
00722 void GsfElectronAlgo::endEvent()
00723  {
00724   if (eventData_==0)
00725    { throw cms::Exception("GsfElectronAlgo|InternalError")<<"lacking event data" ; }
00726   delete eventData_ ;
00727   eventData_ = 0 ;
00728  }
00729 
00730 void GsfElectronAlgo::displayInternalElectrons( const std::string & title ) const
00731  {
00732   LogTrace("GsfElectronAlgo") << "========== " << title << " ==========";
00733   LogTrace("GsfElectronAlgo") << "Event: " << eventData_->event->id();
00734   LogTrace("GsfElectronAlgo") << "Number of electrons: " << eventData_->electrons->size() ;
00735   GsfElectronPtrCollection::const_iterator it ;
00736   for ( it = eventData_->electrons->begin(); it != eventData_->electrons->end(); it++ )
00737    {
00738     LogTrace("GsfElectronAlgo") << "Electron with charge, pt, eta, phi: "  << (*it)->charge() << " , "
00739         << (*it)->pt() << " , " << (*it)->eta() << " , " << (*it)->phi();
00740    }
00741   LogTrace("GsfElectronAlgo") << "=================================================";
00742  }
00743 
00744 void GsfElectronAlgo::completeElectrons()
00745  {
00746   if (electronData_!=0)
00747    { throw cms::Exception("GsfElectronAlgo|InternalError")<<"unexpected electron data" ; }
00748 
00749   const GsfElectronCoreCollection * coreCollection = eventData_->coreElectrons.product() ;
00750   for ( unsigned int i=0 ; i<coreCollection->size() ; ++i )
00751    {
00752     // check there is no existing electron with this core
00753     const GsfElectronCoreRef coreRef = edm::Ref<GsfElectronCoreCollection>(eventData_->coreElectrons,i) ;
00754     bool coreFound = false ;
00755     GsfElectronPtrCollection::const_iterator itrEle ;
00756     for
00757      ( itrEle = eventData_->electrons->begin() ;
00758        itrEle != eventData_->electrons->end() ;
00759        itrEle++ )
00760      {
00761       if ((*itrEle)->core()==coreRef)
00762        {
00763         coreFound = true ;
00764         break ;
00765        }
00766      }
00767     if (coreFound) continue ;
00768 
00769     // check there is a super-cluster
00770     if (coreRef->superCluster().isNull()) continue ;
00771 
00772     // prepare internal structure for electron specific data
00773     delete electronData_ ;
00774     electronData_ = new ElectronData(coreRef,*eventData_->beamspot) ;
00775 
00776     // calculate and check Trajectory StatesOnSurface....
00777     if ( !electronData_->calculateTSOS( eventSetupData_->mtsTransform, eventSetupData_->constraintAtVtx ) ) continue ;
00778 
00779     createElectron() ;
00780 
00781    } // loop over tracks
00782 
00783   delete electronData_ ;
00784   electronData_ = 0 ;
00785  }
00786 
00787 void GsfElectronAlgo::clonePreviousElectrons()
00788  {
00789   const GsfElectronCollection * oldElectrons = eventData_->previousElectrons.product() ;
00790   const GsfElectronCoreCollection * newCores = eventData_->coreElectrons.product() ;
00791   GsfElectronCollection::const_iterator oldElectron ;
00792   for
00793    ( oldElectron = oldElectrons->begin() ;
00794      oldElectron != oldElectrons->end() ;
00795      ++oldElectron )
00796    {
00797     const GsfElectronCoreRef oldCoreRef = oldElectron->core() ;
00798     const GsfTrackRef oldElectronGsfTrackRef = oldCoreRef->gsfTrack() ;
00799     unsigned int icore ;
00800     for ( icore=0 ; icore<newCores->size() ; ++icore )
00801      {
00802       if (oldElectronGsfTrackRef==(*newCores)[icore].gsfTrack())
00803        {
00804         const GsfElectronCoreRef coreRef = edm::Ref<GsfElectronCoreCollection>(eventData_->coreElectrons,icore) ;
00805         eventData_->electrons->push_back(new GsfElectron(*oldElectron,coreRef)) ;
00806         break ;
00807        }
00808      }
00809    }
00810  }
00811 
00812 void GsfElectronAlgo::addPflowInfo()
00813  {
00814   bool found ;
00815   const GsfElectronCollection * pfElectrons = eventData_->pflowElectrons.product() ;
00816   GsfElectronCollection::const_iterator pfElectron ;
00817 
00818   GsfElectronPtrCollection::iterator el ;
00819   for
00820    ( el = eventData_->electrons->begin() ;
00821      el != eventData_->electrons->end() ;
00822      el++ )
00823    {
00824 //    // MVA
00825 //    // we check that the value is never inferior to the no-cut value
00826 //    // we generally use in the configuration file for minMVA.
00827 //    GsfTrackRef gsfTrackRef = (*el)->gsfTrack() ;
00828 //    float mva = (*eventData_->pfMva.product())[gsfTrackRef] ;
00829 //    if (mva<noCutMin) { throw cms::Exception("GsfElectronAlgo|UnexpectedMvaValue")<<"unexpected MVA value: "<<mva ; }
00830 //
00831 //    // Mva Output
00832 //    GsfElectron::MvaOutput mvaOutput ;
00833 //    mvaOutput.mva = mva ;
00834 //    (*el)->setMvaOutput(mvaOutput) ;
00835 
00836     // Retreive info from pflow electrons
00837     found = false ;
00838     for
00839      ( pfElectron = pfElectrons->begin() ; pfElectron != pfElectrons->end() ; pfElectron++ )
00840      {
00841       if (pfElectron->gsfTrack()==(*el)->gsfTrack())
00842        {
00843         if (found)
00844          {
00845           edm::LogWarning("GsfElectronProducer")<<"associated pfGsfElectron already found" ;
00846          }
00847         else
00848          {
00849           found = true ;
00850           (*el)->setPfIsolationVariables(pfElectron->pfIsolationVariables()) ;
00851           (*el)->setMvaInput(pfElectron->mvaInput()) ;
00852           (*el)->setMvaOutput(pfElectron->mvaOutput()) ;
00853           double noCutMin = -999999999 ;
00854           if ((*el)->mva()<noCutMin) { throw cms::Exception("GsfElectronAlgo|UnexpectedMvaValue")<<"unexpected MVA value: "<<(*el)->mva() ; }
00855          }
00856        }
00857      }
00858 
00859     // Preselection
00860     setMvaPreselectionFlag(*el) ;
00861 
00862     // Shower Shape of pflow cluster
00863     if (!((*el)->pflowSuperCluster().isNull()))
00864      {
00865       reco::GsfElectron::ShowerShape pflowShowerShape ;
00866       calculateShowerShape((*el)->pflowSuperCluster(),true,pflowShowerShape) ;
00867       (*el)->setPfShowerShape(pflowShowerShape) ;
00868      }
00869     else if ((*el)->passingMvaPreselection())
00870     { edm::LogError("GsfElectronCoreProducer")<<"Preselected tracker driven GsfTrack with no associated pflow SuperCluster." ; }
00871    }
00872  }
00873 
00874 bool GsfElectronAlgo::isPreselected( GsfElectron * ele )
00875  { return (ele->passingCutBasedPreselection()||ele->passingMvaPreselection()) ; }
00876 
00877 void GsfElectronAlgo::removeNotPreselectedElectrons()
00878  {
00879   GsfElectronPtrCollection::size_type ei = 1, emax = eventData_->electrons->size() ;
00880   GsfElectronPtrCollection::iterator eitr = eventData_->electrons->begin() ;
00881   while (eitr!=eventData_->electrons->end())
00882    {
00883     LogTrace("GsfElectronAlgo")<<"========== removed not preselected "<<ei<<"/"<<emax<<"==========" ;
00884     if (isPreselected(*eitr))
00885      { ++eitr ; ++ei ; }
00886     else
00887      { delete (*eitr) ; eitr = eventData_->electrons->erase(eitr) ; ++ei ; }
00888    }
00889  }
00890 
00891 void GsfElectronAlgo::setCutBasedPreselectionFlag( GsfElectron * ele, const reco::BeamSpot & bs )
00892  {
00893   // default value
00894   ele->setPassCutBasedPreselection(false) ;
00895 
00896   // kind of seeding
00897   bool eg = ele->core()->ecalDrivenSeed() ;
00898   bool pf = ele->core()->trackerDrivenSeed() && !ele->core()->ecalDrivenSeed() ;
00899   if (eg&&pf) { throw cms::Exception("GsfElectronAlgo|BothEcalAndPureTrackerDriven")<<"An electron cannot be both egamma and purely pflow" ; }
00900   if ((!eg)&&(!pf)) { throw cms::Exception("GsfElectronAlgo|NeitherEcalNorPureTrackerDriven")<<"An electron cannot be neither egamma nor purely pflow" ; }
00901   const CutsConfiguration * cfg = (eg?&generalData_->cutsCfg:&generalData_->cutsCfgPflow) ;
00902 
00903   // Et cut
00904   double etaValue = EleRelPoint(ele->superCluster()->position(),bs.position()).eta() ;
00905   double etValue = ele->superCluster()->energy()/cosh(etaValue) ;
00906   LogTrace("GsfElectronAlgo") << "Et : " << etValue ;
00907   if (ele->isEB() && (etValue < cfg->minSCEtBarrel)) return ;
00908   if (ele->isEE() && (etValue < cfg->minSCEtEndcaps)) return ;
00909   LogTrace("GsfElectronAlgo") << "Et criteria are satisfied";
00910 
00911   // E/p cut
00912   double eopValue = ele->eSuperClusterOverP() ;
00913   LogTrace("GsfElectronAlgo") << "E/p : " << eopValue ;
00914   if (ele->isEB() && (eopValue > cfg->maxEOverPBarrel)) return ;
00915   if (ele->isEE() && (eopValue > cfg->maxEOverPEndcaps)) return ;
00916   if (ele->isEB() && (eopValue < cfg->minEOverPBarrel)) return ;
00917   if (ele->isEE() && (eopValue < cfg->minEOverPEndcaps)) return ;
00918   LogTrace("GsfElectronAlgo") << "E/p criteria are satisfied";
00919 
00920   // HoE cuts
00921   LogTrace("GsfElectronAlgo") << "HoE1 : " << ele->hcalDepth1OverEcal() << ", HoE2 : " << ele->hcalDepth2OverEcal();
00922   double had = ele->hcalOverEcal()*ele->superCluster()->energy() ;
00923   const reco::CaloCluster & seedCluster = *(ele->superCluster()->seed()) ;
00924   int detector = seedCluster.hitsAndFractions()[0].first.subdetId() ;
00925   bool HoEveto = false ;
00926   if (detector==EcalBarrel && (had<cfg->maxHBarrel || (had/ele->superCluster()->energy())<cfg->maxHOverEBarrel)) HoEveto=true;
00927   else if (detector==EcalEndcap && (had<cfg->maxHEndcaps || (had/ele->superCluster()->energy())<cfg->maxHOverEEndcaps)) HoEveto=true;
00928   if ( !HoEveto ) return ;
00929   LogTrace("GsfElectronAlgo") << "H/E criteria are satisfied";
00930 
00931   // delta eta criteria
00932   double deta = ele->deltaEtaSuperClusterTrackAtVtx() ;
00933   LogTrace("GsfElectronAlgo") << "delta eta : " << deta ;
00934   if (ele->isEB() && (std::abs(deta) > cfg->maxDeltaEtaBarrel)) return ;
00935   if (ele->isEE() && (std::abs(deta) > cfg->maxDeltaEtaEndcaps)) return ;
00936   LogTrace("GsfElectronAlgo") << "Delta eta criteria are satisfied";
00937 
00938   // delta phi criteria
00939   double dphi = ele->deltaPhiSuperClusterTrackAtVtx();
00940   LogTrace("GsfElectronAlgo") << "delta phi : " << dphi;
00941   if (ele->isEB() && (std::abs(dphi) > cfg->maxDeltaPhiBarrel)) return ;
00942   if (ele->isEE() && (std::abs(dphi) > cfg->maxDeltaPhiEndcaps)) return ;
00943   LogTrace("GsfElectronAlgo") << "Delta phi criteria are satisfied";
00944 
00945   // sigma ieta ieta
00946   LogTrace("GsfElectronAlgo") << "sigma ieta ieta : " << ele->sigmaIetaIeta();
00947   if (ele->isEB() && (ele->sigmaIetaIeta() > cfg->maxSigmaIetaIetaBarrel)) return ;
00948   if (ele->isEE() && (ele->sigmaIetaIeta() > cfg->maxSigmaIetaIetaEndcaps)) return ;
00949   LogTrace("GsfElectronAlgo") << "Sigma ieta ieta criteria are satisfied";
00950 
00951   // fiducial
00952   if (!ele->isEB() && cfg->isBarrel) return ;
00953   if (!ele->isEE() && cfg->isEndcaps) return ;
00954   if (cfg->isFiducial && (ele->isEBEEGap()||ele->isEBEtaGap()||ele->isEBPhiGap()||ele->isEERingGap()||ele->isEEDeeGap())) return ;
00955   LogTrace("GsfElectronAlgo") << "Fiducial flags criteria are satisfied";
00956 
00957   // seed in TEC
00958   edm::RefToBase<TrajectorySeed> seed = ele->gsfTrack()->extra()->seedRef() ;
00959   ElectronSeedRef elseed = seed.castTo<ElectronSeedRef>() ;
00960   if (eg && !generalData_->cutsCfg.seedFromTEC)
00961    {
00962     if (elseed.isNull())
00963      { throw cms::Exception("GsfElectronAlgo|NotElectronSeed")<<"The GsfTrack seed is not an ElectronSeed ?!" ; }
00964     else
00965      { if (elseed->subDet2()==6) return ; }
00966    }
00967 
00968   // transverse impact parameter
00969   if (std::abs(ele->gsfTrack()->dxy(bs.position()))>cfg->maxTIP) return ;
00970   LogTrace("GsfElectronAlgo") << "TIP criterion is satisfied" ;
00971 
00972   LogTrace("GsfElectronAlgo") << "All cut based criteria are satisfied" ;
00973   ele->setPassCutBasedPreselection(true) ;
00974  }
00975 
00976 void GsfElectronAlgo::setMvaPreselectionFlag( GsfElectron * ele )
00977  {
00978   ele->setPassMvaPreselection(false) ;
00979   if (ele->core()->ecalDrivenSeed())
00980    { if (ele->mva()>=generalData_->cutsCfg.minMVA) ele->setPassMvaPreselection(true) ; }
00981   else
00982    { if (ele->mva()>=generalData_->cutsCfgPflow.minMVA) ele->setPassMvaPreselection(true) ; }
00983   if (ele->passingMvaPreselection())
00984    { LogTrace("GsfElectronAlgo") << "Mva criterion is satisfied" ; }
00985  }
00986 
00987 void GsfElectronAlgo::createElectron()
00988  {
00989   // eventually check ctf track
00990   if (generalData_->strategyCfg.ctfTracksCheck)
00991    { electronData_->checkCtfTrack(eventData_->currentCtfTracks) ; }
00992 
00993   // charge ID
00994   int eleCharge ;
00995   GsfElectron::ChargeInfo eleChargeInfo ;
00996   electronData_->computeCharge(eleCharge,eleChargeInfo) ;
00997 
00998   // electron basic cluster
00999   CaloClusterPtr elbcRef = electronData_->getEleBasicCluster(eventSetupData_->mtsTransform) ;
01000 
01001   // Seed cluster
01002   const reco::CaloCluster & seedCluster = *(electronData_->superClusterRef->seed()) ;
01003 
01004   // seed Xtal
01005   // temporary, till CaloCluster->seed() is made available
01006   DetId seedXtalId = seedCluster.hitsAndFractions()[0].first ;
01007 
01008   electronData_->calculateMode(eventSetupData_->mtsMode) ;
01009 
01010 
01011   //====================================================
01012   // Candidate attributes
01013   //====================================================
01014 
01015   Candidate::LorentzVector momentum = electronData_->calculateMomentum() ;
01016 
01017 
01018   //====================================================
01019   // Track-Cluster Matching
01020   //====================================================
01021 
01022   reco::GsfElectron::TrackClusterMatching tcMatching ;
01023   tcMatching.electronCluster = elbcRef ;
01024   tcMatching.eSuperClusterOverP = (electronData_->vtxMom.mag()>0)?(electronData_->superClusterRef->energy()/electronData_->vtxMom.mag()):(-1.) ;
01025   tcMatching.eSeedClusterOverP = (electronData_->vtxMom.mag()>0.)?(seedCluster.energy()/electronData_->vtxMom.mag()):(-1) ;
01026   tcMatching.eSeedClusterOverPout = (electronData_->seedMom.mag()>0.)?(seedCluster.energy()/electronData_->seedMom.mag()):(-1.) ;
01027   tcMatching.eEleClusterOverPout = (electronData_->eleMom.mag()>0.)?(elbcRef->energy()/electronData_->eleMom.mag()):(-1.) ;
01028 
01029   EleRelPointPair scAtVtx(electronData_->superClusterRef->position(),electronData_->sclPos,eventData_->beamspot->position()) ;
01030   tcMatching.deltaEtaSuperClusterAtVtx = scAtVtx.dEta() ;
01031   tcMatching.deltaPhiSuperClusterAtVtx = scAtVtx.dPhi() ;
01032 
01033   EleRelPointPair seedAtCalo(seedCluster.position(),electronData_->seedPos,eventData_->beamspot->position()) ;
01034   tcMatching.deltaEtaSeedClusterAtCalo = seedAtCalo.dEta() ;
01035   tcMatching.deltaPhiSeedClusterAtCalo = seedAtCalo.dPhi() ;
01036 
01037   EleRelPointPair ecAtCalo(elbcRef->position(),electronData_->elePos,eventData_->beamspot->position()) ;
01038   tcMatching.deltaEtaEleClusterAtCalo = ecAtCalo.dEta() ;
01039   tcMatching.deltaPhiEleClusterAtCalo = ecAtCalo.dPhi() ;
01040 
01041 
01042   //=======================================================
01043   // Track extrapolations
01044   //=======================================================
01045 
01046   reco::GsfElectron::TrackExtrapolations tkExtra ;
01047   ele_convert(electronData_->vtxPos,tkExtra.positionAtVtx) ;
01048   ele_convert(electronData_->sclPos,tkExtra.positionAtCalo) ;
01049   ele_convert(electronData_->vtxMom,tkExtra.momentumAtVtx) ;
01050   ele_convert(electronData_->sclMom,tkExtra.momentumAtCalo) ;
01051   ele_convert(electronData_->seedMom,tkExtra.momentumOut) ;
01052   ele_convert(electronData_->eleMom,tkExtra.momentumAtEleClus) ;
01053   ele_convert(electronData_->vtxMomWithConstraint,tkExtra.momentumAtVtxWithConstraint) ;
01054 
01055 
01056   //=======================================================
01057   // Closest Ctf Track
01058   //=======================================================
01059 
01060   reco::GsfElectron::ClosestCtfTrack ctfInfo ;
01061   ctfInfo.ctfTrack = electronData_->ctfTrackRef  ;
01062   ctfInfo.shFracInnerHits = electronData_->shFracInnerHits ;
01063 
01064 
01065   //====================================================
01066   // FiducialFlags, using nextToBoundary definition of gaps
01067   //====================================================
01068 
01069   reco::GsfElectron::FiducialFlags fiducialFlags ;
01070   int detector = seedXtalId.subdetId() ;
01071   double feta=std::abs(electronData_->superClusterRef->position().eta()) ;
01072   if (detector==EcalBarrel)
01073    {
01074     fiducialFlags.isEB = true ;
01075     EBDetId ebdetid(seedXtalId);
01076     if (EBDetId::isNextToEtaBoundary(ebdetid))
01077      {
01078       if (ebdetid.ietaAbs()==85)
01079        { fiducialFlags.isEBEEGap = true ; }
01080       else
01081        { fiducialFlags.isEBEtaGap = true ; }
01082      }
01083     if (EBDetId::isNextToPhiBoundary(ebdetid))
01084      { fiducialFlags.isEBPhiGap = true ; }
01085    }
01086   else if (detector==EcalEndcap)
01087    {
01088     fiducialFlags.isEE = true ;
01089     EEDetId eedetid(seedXtalId);
01090     if (EEDetId::isNextToRingBoundary(eedetid))
01091      {
01092       if (std::abs(feta)<2.)
01093        { fiducialFlags.isEBEEGap = true ; }
01094       else
01095        { fiducialFlags.isEERingGap = true ; }
01096      }
01097     if (EEDetId::isNextToDBoundary(eedetid))
01098      { fiducialFlags.isEEDeeGap = true ; }
01099    }
01100   else
01101    { throw cms::Exception("GsfElectronAlgo|UnknownXtalRegion")<<"createElectron(): do not know if it is a barrel or endcap seed cluster !!!!" ; }
01102 
01103 
01104   //====================================================
01105   // ShowerShape
01106   //====================================================
01107 
01108   reco::GsfElectron::ShowerShape showerShape ;
01109   calculateShowerShape(electronData_->superClusterRef,!(electronData_->coreRef->ecalDrivenSeed()),showerShape) ;
01110 
01111 
01112   //====================================================
01113   // ConversionRejection
01114   //====================================================
01115 
01116   eventData_->retreiveOriginalTrackCollections(electronData_->ctfTrackRef,electronData_->coreRef->gsfTrack()) ;
01117 
01118   ConversionFinder conversionFinder ;
01119   double BInTesla = eventSetupData_->magField->inTesla(GlobalPoint(0.,0.,0.)).z() ;
01120   edm::Handle<reco::TrackCollection> ctfTracks = eventData_->originalCtfTracks ;
01121   if (!ctfTracks.isValid()) { ctfTracks = eventData_->currentCtfTracks ; }
01122 
01123   // values of conversionInfo.flag()
01124   // -9999 : Partner track was not found
01125   // 0     : Partner track found in the CTF collection using
01126   // 1     : Partner track found in the CTF collection using
01127   // 2     : Partner track found in the GSF collection using
01128   // 3     : Partner track found in the GSF collection using the electron's GSF track
01129   ConversionInfo conversionInfo = conversionFinder.getConversionInfo
01130    (*electronData_->coreRef,ctfTracks,eventData_->originalGsfTracks,BInTesla) ;
01131 
01132   reco::GsfElectron::ConversionRejection conversionVars ;
01133   conversionVars.flags = conversionInfo.flag()  ;
01134   conversionVars.dist = conversionInfo.dist()  ;
01135   conversionVars.dcot = conversionInfo.dcot()  ;
01136   conversionVars.radius = conversionInfo.radiusOfConversion()  ;
01137   if ((conversionVars.flags==0)or(conversionVars.flags==1))
01138     conversionVars.partner = TrackBaseRef(conversionInfo.conversionPartnerCtfTk())  ;
01139   else if ((conversionVars.flags==2)or(conversionVars.flags==3))
01140     conversionVars.partner = TrackBaseRef(conversionInfo.conversionPartnerGsfTk())  ;
01141 
01142 
01143   //====================================================
01144   // brems fraction
01145   //====================================================
01146 
01147   float fbrem = (electronData_->outMom.mag()>0.)?((electronData_->innMom.mag()-electronData_->outMom.mag())/electronData_->innMom.mag()):1.e30 ;
01148 
01149 
01150   //====================================================
01151   // Go !
01152   //====================================================
01153 
01154   GsfElectron * ele = new
01155     GsfElectron
01156      ( eleCharge,eleChargeInfo,electronData_->coreRef,
01157        tcMatching, tkExtra, ctfInfo,
01158        fiducialFlags,showerShape,
01159        conversionVars,
01160        fbrem) ;
01161   ele->setP4(GsfElectron::P4_FROM_SUPER_CLUSTER,momentum,0,true) ;
01162 
01163   // set corrections + classification
01164   ElectronClassification theClassifier;
01165   theClassifier.correct(*ele);
01166   // energy corrections only for ecalDriven electrons
01167   if (ele->core()->ecalDrivenSeed()) {
01168     ElectronEnergyCorrector theEnCorrector(generalData_->superClusterErrorFunction);
01169     theEnCorrector.correct(*ele, *eventData_->beamspot, generalData_->strategyCfg.applyEtaCorrection);
01170     ElectronMomentumCorrector theMomCorrector;
01171     theMomCorrector.correct(*ele,electronData_->vtxTSOS);
01172   }
01173 
01174 
01175   // now isolation variables
01176   reco::GsfElectron::IsolationVariables dr03, dr04 ;
01177   dr03.tkSumPt = eventData_->tkIsolation03->getPtTracks(ele);
01178   dr03.hcalDepth1TowerSumEt = eventData_->hadDepth1Isolation03->getTowerEtSum(ele);
01179   dr03.hcalDepth2TowerSumEt = eventData_->hadDepth2Isolation03->getTowerEtSum(ele);
01180   dr03.ecalRecHitSumEt = eventData_->ecalBarrelIsol03->getEtSum(ele)+eventData_->ecalEndcapIsol03->getEtSum(ele);
01181   dr04.tkSumPt = eventData_->tkIsolation04->getPtTracks(ele);
01182   dr04.hcalDepth1TowerSumEt = eventData_->hadDepth1Isolation04->getTowerEtSum(ele);
01183   dr04.hcalDepth2TowerSumEt = eventData_->hadDepth2Isolation04->getTowerEtSum(ele);
01184   dr04.ecalRecHitSumEt = eventData_->ecalBarrelIsol04->getEtSum(ele)+eventData_->ecalEndcapIsol04->getEtSum(ele);
01185   ele->setIsolation03(dr03);
01186   ele->setIsolation04(dr04);
01187 
01188   // preselection flag
01189   setCutBasedPreselectionFlag(ele,*eventData_->beamspot) ;
01190 
01191   LogTrace("GsfElectronAlgo")<<"Constructed new electron with energy  "<< ele->p4().e() ;
01192 
01193   eventData_->electrons->push_back(ele) ;
01194  }
01195 
01196 
01197 //=======================================================================================
01198 // Ambiguity solving
01199 //=======================================================================================
01200 
01201 //bool better_electron( const reco::GsfElectron * e1, const reco::GsfElectron * e2 )
01202 // { return (std::abs(e1->eSuperClusterOverP()-1)<std::abs(e2->eSuperClusterOverP()-1)) ; }
01203 
01204 void GsfElectronAlgo::setAmbiguityData( bool ignoreNotPreselected )
01205  {
01206   GsfElectronPtrCollection::iterator e1, e2 ;
01207   if (generalData_->strategyCfg.ambSortingStrategy==0)
01208    { eventData_->electrons->sort(EgAmbiguityTools::isBetter) ; }
01209   else if (generalData_->strategyCfg.ambSortingStrategy==1)
01210    { eventData_->electrons->sort(EgAmbiguityTools::isInnerMost(eventSetupData_->trackerHandle)) ; }
01211   else
01212    { throw cms::Exception("GsfElectronAlgo|UnknownAmbiguitySortingStrategy")<<"value of generalData_->strategyCfg.ambSortingStrategy is : "<<generalData_->strategyCfg.ambSortingStrategy ; }
01213 
01214   // init
01215   for
01216    ( e1 = eventData_->electrons->begin() ;
01217      e1 != eventData_->electrons->end() ;
01218      ++e1 )
01219    {
01220     (*e1)->clearAmbiguousGsfTracks() ;
01221     (*e1)->setAmbiguous(false) ;
01222    }
01223 
01224   // get ambiguous from GsfPfRecTracks
01225   if (generalData_->strategyCfg.useGsfPfRecTracks)
01226    {
01227     for
01228      ( e1 = eventData_->electrons->begin() ;
01229        e1 != eventData_->electrons->end() ;
01230        ++e1 )
01231      {
01232       bool found = false ;
01233       const GsfPFRecTrackCollection * gsfPfRecTrackCollection = eventData_->gsfPfRecTracks.product() ;
01234       GsfPFRecTrackCollection::const_iterator gsfPfRecTrack ;
01235       for ( gsfPfRecTrack=gsfPfRecTrackCollection->begin() ;
01236             gsfPfRecTrack!=gsfPfRecTrackCollection->end() ;
01237             ++gsfPfRecTrack )
01238        {
01239         if (gsfPfRecTrack->gsfTrackRef()==(*e1)->gsfTrack())
01240          {
01241           if (found)
01242            {
01243             edm::LogWarning("GsfElectronAlgo")<<"associated gsfPfRecTrack already found" ;
01244            }
01245           else
01246            {
01247             found = true ;
01248             const std::vector<reco::GsfPFRecTrackRef> & duplicates(gsfPfRecTrack->convBremGsfPFRecTrackRef()) ;
01249             std::vector<reco::GsfPFRecTrackRef>::const_iterator duplicate ;
01250             for ( duplicate = duplicates.begin() ; duplicate != duplicates.end() ; duplicate ++ )
01251              { (*e1)->addAmbiguousGsfTrack((*duplicate)->gsfTrackRef()) ; }
01252            }
01253          }
01254        }
01255      }
01256    }
01257   // or search overlapping clusters
01258   else
01259    {
01260     for
01261      ( e1 = eventData_->electrons->begin() ;
01262        e1 != eventData_->electrons->end() ;
01263        ++e1 )
01264      {
01265       if ((*e1)->ambiguous()) continue ;
01266       if ( ignoreNotPreselected && !isPreselected(*e1) ) continue ;
01267 
01268       SuperClusterRef scRef1 = (*e1)->superCluster();
01269       CaloClusterPtr eleClu1 = (*e1)->electronCluster();
01270       LogDebug("GsfElectronAlgo")
01271         << "Blessing electron with E/P " << (*e1)->eSuperClusterOverP()
01272         << ", cluster " << scRef1.get()
01273         << " & track " << (*e1)->gsfTrack().get() ;
01274 
01275       for
01276        ( e2 = e1, ++e2 ;
01277          e2 != eventData_->electrons->end() ;
01278          ++e2 )
01279        {
01280         if ((*e2)->ambiguous()) continue ;
01281         if ( ignoreNotPreselected && !isPreselected(*e2) ) continue ;
01282 
01283         SuperClusterRef scRef2 = (*e2)->superCluster();
01284         CaloClusterPtr eleClu2 = (*e2)->electronCluster();
01285 
01286         // search if same cluster
01287         bool sameCluster = false ;
01288         if (generalData_->strategyCfg.ambClustersOverlapStrategy==0)
01289          { sameCluster = (scRef1==scRef2) ; }
01290         else if (generalData_->strategyCfg.ambClustersOverlapStrategy==1)
01291          {
01292           float eMin = 1. ;
01293           float threshold = eMin*cosh(EleRelPoint(scRef1->position(),eventData_->beamspot->position()).eta()) ;
01294           sameCluster =
01295            ( (EgAmbiguityTools::sharedEnergy(&(*eleClu1),&(*eleClu2),eventData_->reducedEBRecHits,eventData_->reducedEERecHits)>=threshold) ||
01296              (EgAmbiguityTools::sharedEnergy(&(*scRef1->seed()),&(*eleClu2),eventData_->reducedEBRecHits,eventData_->reducedEERecHits)>=threshold) ||
01297              (EgAmbiguityTools::sharedEnergy(&(*eleClu1),&(*scRef2->seed()),eventData_->reducedEBRecHits,eventData_->reducedEERecHits)>=threshold) ||
01298              (EgAmbiguityTools::sharedEnergy(&(*scRef1->seed()),&(*scRef2->seed()),eventData_->reducedEBRecHits,eventData_->reducedEERecHits)>=threshold) ) ;
01299          }
01300         else
01301          { throw cms::Exception("GsfElectronAlgo|UnknownAmbiguityClustersOverlapStrategy")<<"value of generalData_->strategyCfg.ambClustersOverlapStrategy is : "<<generalData_->strategyCfg.ambClustersOverlapStrategy ; }
01302 
01303         // main instructions
01304         if (sameCluster)
01305          {
01306           LogDebug("GsfElectronAlgo")
01307             << "Discarding electron with E/P " << (*e2)->eSuperClusterOverP()
01308             << ", cluster " << scRef2.get()
01309             << " and track " << (*e2)->gsfTrack().get() ;
01310           (*e1)->addAmbiguousGsfTrack((*e2)->gsfTrack()) ;
01311           (*e2)->setAmbiguous(true) ;
01312          }
01313         else if ((*e1)->gsfTrack()==(*e2)->gsfTrack())
01314          {
01315           edm::LogWarning("GsfElectronAlgo")
01316             << "Forgetting electron with E/P " << (*e2)->eSuperClusterOverP()
01317             << ", cluster " << scRef2.get()
01318             << " and track " << (*e2)->gsfTrack().get() ;
01319           (*e2)->setAmbiguous(true) ;
01320          }
01321        }
01322      }
01323    }
01324  }
01325 
01326 void GsfElectronAlgo::removeAmbiguousElectrons()
01327  {
01328   GsfElectronPtrCollection::size_type ei = 1, emax = eventData_->electrons->size() ;
01329   GsfElectronPtrCollection::iterator eitr = eventData_->electrons->begin() ;
01330   while (eitr!=eventData_->electrons->end())
01331    {
01332     LogTrace("GsfElectronAlgo")<<"========== remove ambiguous "<<ei<<"/"<<emax<<"==========" ;
01333     if ((*eitr)->ambiguous())
01334      { delete (*eitr) ; eitr = eventData_->electrons->erase(eitr) ; ++ei ; }
01335     else
01336      { ++eitr ; ++ei ; }
01337    }
01338  }
01339
/data/doxygen/doxygen-1.7.3/gen/CMSSW_4_2_8/src/RecoEgamma/EgammaElectronAlgos/src/GsfElectronAlgo.cc
