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