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