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