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