Reference Manual

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