CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_6_1_2_SLHC4_patch1/src/RecoEgamma/EgammaElectronAlgos/src/GsfElectronAlgo.cc

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