00001
00002
00003
00004
00005 #include "RecoParticleFlow/PFProducer/interface/PFElectronAlgo.h"
00006 #include "RecoParticleFlow/PFProducer/interface/PFMuonAlgo.h"
00007
00008 #include "DataFormats/ParticleFlowReco/interface/PFBlockElementCluster.h"
00009 #include "DataFormats/ParticleFlowReco/interface/PFBlockElementBrem.h"
00010 #include "DataFormats/ParticleFlowReco/interface/PFRecHitFraction.h"
00011 #include "DataFormats/ParticleFlowReco/interface/PFRecHitFwd.h"
00012 #include "DataFormats/ParticleFlowReco/interface/PFRecHit.h"
00013 #include "DataFormats/ParticleFlowReco/interface/PFCluster.h"
00014 #include "DataFormats/GsfTrackReco/interface/GsfTrack.h"
00015 #include "DataFormats/GsfTrackReco/interface/GsfTrackFwd.h"
00016 #include "RecoParticleFlow/PFProducer/interface/PFElectronExtraEqual.h"
00017 #include "RecoParticleFlow/PFClusterTools/interface/PFEnergyCalibration.h"
00018 #include "RecoParticleFlow/PFClusterTools/interface/PFSCEnergyCalibration.h"
00019 #include "RecoParticleFlow/PFClusterTools/interface/PFEnergyResolution.h"
00020 #include "RecoParticleFlow/PFClusterTools/interface/PFClusterWidthAlgo.h"
00021 #include "DataFormats/ParticleFlowReco/interface/PFRecTrack.h"
00022 #include "DataFormats/ParticleFlowReco/interface/GsfPFRecTrack.h"
00023
00024 #include <iomanip>
00025 #include <algorithm>
00026
00027 using namespace std;
00028 using namespace reco;
00029 PFElectronAlgo::PFElectronAlgo(const double mvaEleCut,
00030 string mvaWeightFileEleID,
00031 const boost::shared_ptr<PFSCEnergyCalibration>& thePFSCEnergyCalibration,
00032 bool applyCrackCorrections,
00033 bool usePFSCEleCalib,
00034 bool useEGElectrons,
00035 bool useEGammaSupercluster,
00036 double sumEtEcalIsoForEgammaSC_barrel,
00037 double sumEtEcalIsoForEgammaSC_endcap,
00038 double coneEcalIsoForEgammaSC,
00039 double sumPtTrackIsoForEgammaSC_barrel,
00040 double sumPtTrackIsoForEgammaSC_endcap,
00041 unsigned int nTrackIsoForEgammaSC,
00042 double coneTrackIsoForEgammaSC):
00043 mvaEleCut_(mvaEleCut),
00044 thePFSCEnergyCalibration_(thePFSCEnergyCalibration),
00045 applyCrackCorrections_(applyCrackCorrections),
00046 usePFSCEleCalib_(usePFSCEleCalib),
00047 useEGElectrons_(useEGElectrons),
00048 useEGammaSupercluster_(useEGammaSupercluster),
00049 sumEtEcalIsoForEgammaSC_barrel_(sumEtEcalIsoForEgammaSC_barrel),
00050 sumEtEcalIsoForEgammaSC_endcap_(sumEtEcalIsoForEgammaSC_endcap),
00051 coneEcalIsoForEgammaSC_(coneEcalIsoForEgammaSC),
00052 sumPtTrackIsoForEgammaSC_barrel_(sumPtTrackIsoForEgammaSC_barrel),
00053 sumPtTrackIsoForEgammaSC_endcap_(sumPtTrackIsoForEgammaSC_endcap),
00054 nTrackIsoForEgammaSC_(nTrackIsoForEgammaSC),
00055 coneTrackIsoForEgammaSC_(coneTrackIsoForEgammaSC)
00056 {
00057
00058 tmvaReader_ = new TMVA::Reader("!Color:Silent");
00059 tmvaReader_->AddVariable("lnPt_gsf",&lnPt_gsf);
00060 tmvaReader_->AddVariable("Eta_gsf",&Eta_gsf);
00061 tmvaReader_->AddVariable("dPtOverPt_gsf",&dPtOverPt_gsf);
00062 tmvaReader_->AddVariable("DPtOverPt_gsf",&DPtOverPt_gsf);
00063
00064 tmvaReader_->AddVariable("chi2_gsf",&chi2_gsf);
00065
00066 tmvaReader_->AddVariable("nhit_kf",&nhit_kf);
00067 tmvaReader_->AddVariable("chi2_kf",&chi2_kf);
00068 tmvaReader_->AddVariable("EtotPinMode",&EtotPinMode);
00069 tmvaReader_->AddVariable("EGsfPoutMode",&EGsfPoutMode);
00070 tmvaReader_->AddVariable("EtotBremPinPoutMode",&EtotBremPinPoutMode);
00071 tmvaReader_->AddVariable("DEtaGsfEcalClust",&DEtaGsfEcalClust);
00072 tmvaReader_->AddVariable("SigmaEtaEta",&SigmaEtaEta);
00073 tmvaReader_->AddVariable("HOverHE",&HOverHE);
00074
00075 tmvaReader_->AddVariable("lateBrem",&lateBrem);
00076 tmvaReader_->AddVariable("firstBrem",&firstBrem);
00077 tmvaReader_->BookMVA("BDT",mvaWeightFileEleID.c_str());
00078 }
00079 void PFElectronAlgo::RunPFElectron(const reco::PFBlockRef& blockRef,
00080 std::vector<bool>& active){
00081
00082
00083 AssMap associatedToGsf;
00084 AssMap associatedToBrems;
00085 AssMap associatedToEcal;
00086
00087
00088 elCandidate_.clear();
00089 electronExtra_.clear();
00090 allElCandidate_.clear();
00091 electronConstituents_.clear();
00092 fifthStepKfTrack_.clear();
00093 convGsfTrack_.clear();
00094
00095
00096 bool blockHasGSF = SetLinks(blockRef,associatedToGsf,
00097 associatedToBrems,associatedToEcal,
00098 active);
00099
00100
00101 if (blockHasGSF) {
00102
00103 BDToutput_.clear();
00104
00105 lockExtraKf_.clear();
00106
00107 BDToutput_.assign(associatedToGsf.size(),-1.);
00108 lockExtraKf_.assign(associatedToGsf.size(),true);
00109
00110
00111 SetIDOutputs(blockRef,associatedToGsf,associatedToBrems,associatedToEcal);
00112
00113
00114
00115
00116 SetCandidates(blockRef,associatedToGsf,associatedToBrems,associatedToEcal);
00117 if (elCandidate_.size() > 0 ){
00118 isvalid_ = true;
00119
00120
00121 SetActive(blockRef,associatedToGsf,associatedToBrems,associatedToEcal,active);
00122 }
00123 }
00124 }
00125 bool PFElectronAlgo::SetLinks(const reco::PFBlockRef& blockRef,
00126 AssMap& associatedToGsf_,
00127 AssMap& associatedToBrems_,
00128 AssMap& associatedToEcal_,
00129 std::vector<bool>& active){
00130 unsigned int CutIndex = 100000;
00131 double CutGSFECAL = 10000. ;
00132
00133 PFEnergyCalibration pfcalib_;
00134 bool DebugSetLinksSummary = false;
00135 bool DebugSetLinksDetailed = false;
00136
00137 const reco::PFBlock& block = *blockRef;
00138 const edm::OwnVector< reco::PFBlockElement >& elements = block.elements();
00139 PFBlock::LinkData linkData = block.linkData();
00140
00141 bool IsThereAGSFTrack = false;
00142 bool IsThereAGoodGSFTrack = false;
00143
00144 vector<unsigned int> trackIs(0);
00145 vector<unsigned int> gsfIs(0);
00146 vector<unsigned int> ecalIs(0);
00147
00148 std::vector<bool> localactive(elements.size(),true);
00149
00150
00151
00152 std::multimap<double, unsigned int> kfElems;
00153 for(unsigned int iEle=0; iEle<elements.size(); iEle++) {
00154 localactive[iEle] = active[iEle];
00155 bool thisIsAMuon = false;
00156 PFBlockElement::Type type = elements[iEle].type();
00157 switch( type ) {
00158 case PFBlockElement::TRACK:
00159
00160 thisIsAMuon = PFMuonAlgo::isMuon(elements[iEle]);
00161
00162 if ( !thisIsAMuon && active[iEle] ) {
00163 trackIs.push_back( iEle );
00164 if (DebugSetLinksDetailed)
00165 cout<<"TRACK, stored index, continue "<< iEle << endl;
00166 }
00167 continue;
00168 case PFBlockElement::GSF:
00169
00170 block.associatedElements( iEle,linkData,
00171 kfElems,
00172 reco::PFBlockElement::TRACK,
00173 reco::PFBlock::LINKTEST_ALL );
00174 thisIsAMuon = kfElems.size() ?
00175 PFMuonAlgo::isMuon(elements[kfElems.begin()->second]) : false;
00176
00177 if ( !thisIsAMuon && active[iEle] ) {
00178 IsThereAGSFTrack = true;
00179 gsfIs.push_back( iEle );
00180 if (DebugSetLinksDetailed)
00181 cout<<"GSF, stored index, continue "<< iEle << endl;
00182 }
00183 continue;
00184 case PFBlockElement::ECAL:
00185 if ( active[iEle] ) {
00186 ecalIs.push_back( iEle );
00187 if (DebugSetLinksDetailed)
00188 cout<<"ECAL, stored index, continue "<< iEle << endl;
00189 }
00190 continue;
00191 default:
00192 continue;
00193 }
00194 }
00195
00196
00197 if(IsThereAGSFTrack) {
00198
00199
00200
00201
00202
00203
00204
00205
00206
00207
00208
00209
00210
00211
00212
00213
00214 for(unsigned int iEle=0; iEle<trackIs.size(); iEle++) {
00215 std::multimap<double, unsigned int> gsfElems;
00216 block.associatedElements( trackIs[iEle], linkData,
00217 gsfElems ,
00218 reco::PFBlockElement::GSF,
00219 reco::PFBlock::LINKTEST_ALL );
00220 if(gsfElems.size() == 0){
00221
00222
00223 std::multimap<double, unsigned int> ecalKfElems;
00224 block.associatedElements( trackIs[iEle],linkData,
00225 ecalKfElems,
00226 reco::PFBlockElement::ECAL,
00227 reco::PFBlock::LINKTEST_ALL );
00228 if(ecalKfElems.size() > 0) {
00229 unsigned int ecalKf_index = ecalKfElems.begin()->second;
00230 if(localactive[ecalKf_index]==true) {
00231
00232
00233
00234 bool isGsfLinked = false;
00235 for(unsigned int iGsf=0; iGsf<gsfIs.size(); iGsf++) {
00236
00237
00238
00239 const reco::PFBlockElementGsfTrack * GsfEl =
00240 dynamic_cast<const reco::PFBlockElementGsfTrack*>((&elements[gsfIs[iGsf]]));
00241 if(GsfEl->trackType(reco::PFBlockElement::T_FROM_GAMMACONV)) continue;
00242
00243 std::multimap<double, unsigned int> ecalGsfElems;
00244 block.associatedElements( gsfIs[iGsf],linkData,
00245 ecalGsfElems,
00246 reco::PFBlockElement::ECAL,
00247 reco::PFBlock::LINKTEST_ALL );
00248 if(ecalGsfElems.size() > 0) {
00249 if (ecalGsfElems.begin()->second == ecalKf_index) {
00250 isGsfLinked = true;
00251 }
00252 }
00253 }
00254 if(isGsfLinked == false) {
00255
00256
00257 const reco::PFBlockElementTrack * kfEle =
00258 dynamic_cast<const reco::PFBlockElementTrack*>((&elements[(trackIs[iEle])]));
00259 reco::TrackRef refKf = kfEle->trackRef();
00260 unsigned int Algo = 0;
00261 if (refKf.isNonnull())
00262 Algo = refKf->algo();
00263 if(Algo < 9) {
00264 localactive[ecalKf_index] = false;
00265 }
00266 else {
00267 fifthStepKfTrack_.push_back(make_pair(ecalKf_index,trackIs[iEle]));
00268 }
00269 }
00270 }
00271 }
00272 }
00273 }
00274
00275
00276
00277 for(unsigned int iEle=0; iEle<gsfIs.size(); iEle++) {
00278
00279 if (!localactive[(gsfIs[iEle])]) continue;
00280
00281 localactive[gsfIs[iEle]] = false;
00282 bool ClosestEcalWithKf = false;
00283
00284 if (DebugSetLinksDetailed) cout << " Gsf Index " << gsfIs[iEle] << endl;
00285
00286 const reco::PFBlockElementGsfTrack * GsfEl =
00287 dynamic_cast<const reco::PFBlockElementGsfTrack*>((&elements[(gsfIs[iEle])]));
00288
00289
00290 if(GsfEl->trackType(reco::PFBlockElement::T_FROM_GAMMACONV)) continue;
00291 IsThereAGoodGSFTrack = true;
00292 float eta_gsf = GsfEl->positionAtECALEntrance().eta();
00293 float etaOut_gsf = GsfEl->Pout().eta();
00294 float diffOutEcalEta = fabs(eta_gsf-etaOut_gsf);
00295 reco::GsfTrackRef RefGSF = GsfEl->GsftrackRef();
00296 float Pin_gsf = 0.01;
00297 if (RefGSF.isNonnull() )
00298 Pin_gsf = RefGSF->pMode();
00299
00300
00301
00302 unsigned int KfGsf_index = CutIndex;
00303 unsigned int KfGsf_secondIndex = CutIndex;
00304 std::multimap<double, unsigned int> kfElems;
00305 block.associatedElements( gsfIs[iEle],linkData,
00306 kfElems,
00307 reco::PFBlockElement::TRACK,
00308 reco::PFBlock::LINKTEST_ALL );
00309 std::multimap<double, unsigned int> ecalKfElems;
00310 if (kfElems.size() > 0) {
00311
00312
00313
00314 for(std::multimap<double, unsigned int>::iterator itkf = kfElems.begin();
00315 itkf != kfElems.end(); ++itkf) {
00316 const reco::PFBlockElementTrack * TrkEl =
00317 dynamic_cast<const reco::PFBlockElementTrack*>((&elements[itkf->second]));
00318
00319 bool isPrim = isPrimaryTrack(*TrkEl,*GsfEl);
00320 if(!isPrim)
00321 continue;
00322
00323 if(localactive[itkf->second] == true) {
00324
00325 KfGsf_index = itkf->second;
00326 localactive[KfGsf_index] = false;
00327
00328 block.associatedElements( KfGsf_index, linkData,
00329 ecalKfElems ,
00330 reco::PFBlockElement::ECAL,
00331 reco::PFBlock::LINKTEST_ALL );
00332 }
00333 else {
00334 KfGsf_secondIndex = itkf->second;
00335 }
00336 }
00337 }
00338
00339
00340 std::multimap<double, unsigned int> ecalGsfElems;
00341 block.associatedElements( gsfIs[iEle],linkData,
00342 ecalGsfElems,
00343 reco::PFBlockElement::ECAL,
00344 reco::PFBlock::LINKTEST_ALL );
00345 double ecalGsf_dist = CutGSFECAL;
00346 unsigned int ClosestEcalGsf_index = CutIndex;
00347 if (ecalGsfElems.size() > 0) {
00348 if(localactive[(ecalGsfElems.begin()->second)] == true) {
00349
00350 bool compatibleEPout = true;
00351 if(diffOutEcalEta > 0.3) {
00352 reco::PFClusterRef clusterRef = elements[(ecalGsfElems.begin()->second)].clusterRef();
00353 float EoPout = (clusterRef->energy())/(GsfEl->Pout().t());
00354 if(EoPout > 5)
00355 compatibleEPout = false;
00356 }
00357 if(compatibleEPout) {
00358 ClosestEcalGsf_index = ecalGsfElems.begin()->second;
00359 ecalGsf_dist = block.dist(gsfIs[iEle],ClosestEcalGsf_index,
00360 linkData,reco::PFBlock::LINKTEST_ALL);
00361
00362
00363
00364 std::multimap<double, unsigned int> ecalOtherGsfElems;
00365 block.associatedElements( ClosestEcalGsf_index,linkData,
00366 ecalOtherGsfElems,
00367 reco::PFBlockElement::GSF,
00368 reco::PFBlock::LINKTEST_ALL);
00369
00370 if(ecalOtherGsfElems.size()>0) {
00371
00372 const reco::PFBlockElementGsfTrack * gsfCheck =
00373 dynamic_cast<const reco::PFBlockElementGsfTrack*>((&elements[ecalOtherGsfElems.begin()->second]));
00374
00375 if(ecalOtherGsfElems.begin()->second != gsfIs[iEle]&&
00376 gsfCheck->trackType(reco::PFBlockElement::T_FROM_GAMMACONV) == false) {
00377 ecalGsf_dist = CutGSFECAL;
00378 ClosestEcalGsf_index = CutIndex;
00379 }
00380 }
00381 }
00382
00383 }
00384 }
00385
00386 else if(ecalKfElems.size() > 0) {
00387 if(localactive[(ecalKfElems.begin()->second)] == true) {
00388 ClosestEcalGsf_index = ecalKfElems.begin()->second;
00389 ecalGsf_dist = block.dist(gsfIs[iEle],ClosestEcalGsf_index,
00390 linkData,reco::PFBlock::LINKTEST_ALL);
00391 ClosestEcalWithKf = true;
00392
00393
00394 std::multimap<double, unsigned int> ecalOtherGsfElems;
00395 block.associatedElements( ClosestEcalGsf_index,linkData,
00396 ecalOtherGsfElems,
00397 reco::PFBlockElement::GSF,
00398 reco::PFBlock::LINKTEST_ALL);
00399 if(ecalOtherGsfElems.size() > 0) {
00400 const reco::PFBlockElementGsfTrack * gsfCheck =
00401 dynamic_cast<const reco::PFBlockElementGsfTrack*>((&elements[ecalOtherGsfElems.begin()->second]));
00402
00403 if(ecalOtherGsfElems.begin()->second != gsfIs[iEle] &&
00404 gsfCheck->trackType(reco::PFBlockElement::T_FROM_GAMMACONV) == false) {
00405 ecalGsf_dist = CutGSFECAL;
00406 ClosestEcalGsf_index = CutIndex;
00407 ClosestEcalWithKf = false;
00408 }
00409 }
00410 }
00411 }
00412
00413 if (DebugSetLinksDetailed)
00414 cout << " Closest Ecal to the Gsf/Kf: index " << ClosestEcalGsf_index
00415 << " dist " << ecalGsf_dist << endl;
00416
00417
00418
00419
00420 std::multimap<double, unsigned int> bremElems;
00421 block.associatedElements( gsfIs[iEle],linkData,
00422 bremElems,
00423 reco::PFBlockElement::BREM,
00424 reco::PFBlock::LINKTEST_ALL );
00425
00426
00427 multimap<unsigned int,unsigned int> cleanedEcalBremElems;
00428 vector<unsigned int> keyBremIndex(0);
00429 unsigned int latestBrem_trajP = 0;
00430 unsigned int latestBrem_index = CutIndex;
00431 for(std::multimap<double, unsigned int>::iterator ieb = bremElems.begin();
00432 ieb != bremElems.end(); ++ieb ) {
00433 unsigned int brem_index = ieb->second;
00434 if(localactive[brem_index] == false) continue;
00435
00436
00437
00438 std::multimap<double, unsigned int> ecalBremsElems;
00439
00440 block.associatedElements( brem_index, linkData,
00441 ecalBremsElems,
00442 reco::PFBlockElement::ECAL,
00443 reco::PFBlock::LINKTEST_ALL );
00444
00445 for (std::multimap<double, unsigned int>::iterator ie = ecalBremsElems.begin();
00446 ie != ecalBremsElems.end();ie++) {
00447 unsigned int ecalBrem_index = ie->second;
00448 if(localactive[ecalBrem_index] == false) continue;
00449
00450
00451 float ecalBrem_dist = block.dist(brem_index,ecalBrem_index,
00452 linkData,reco::PFBlock::LINKTEST_ALL);
00453
00454
00455 if (ecalBrem_index == ClosestEcalGsf_index && (ecalBrem_dist + 0.0012) > ecalGsf_dist) continue;
00456
00457
00458 std::multimap<double, unsigned int> sortedBremElems;
00459 block.associatedElements( ecalBrem_index,linkData,
00460 sortedBremElems,
00461 reco::PFBlockElement::BREM,
00462 reco::PFBlock::LINKTEST_ALL);
00463
00464 bool isGoodBrem = false;
00465 unsigned int sortedBrem_index = CutIndex;
00466 for (std::multimap<double, unsigned int>::iterator ibs = sortedBremElems.begin();
00467 ibs != sortedBremElems.end();ibs++) {
00468 unsigned int temp_sortedBrem_index = ibs->second;
00469 std::multimap<double, unsigned int> sortedGsfElems;
00470 block.associatedElements( temp_sortedBrem_index,linkData,
00471 sortedGsfElems,
00472 reco::PFBlockElement::GSF,
00473 reco::PFBlock::LINKTEST_ALL);
00474 bool enteredInPrimaryGsf = false;
00475 for (std::multimap<double, unsigned int>::iterator igs = sortedGsfElems.begin();
00476 igs != sortedGsfElems.end();igs++) {
00477 const reco::PFBlockElementGsfTrack * gsfCheck =
00478 dynamic_cast<const reco::PFBlockElementGsfTrack*>((&elements[igs->second]));
00479
00480 if(gsfCheck->trackType(reco::PFBlockElement::T_FROM_GAMMACONV) == false) {
00481 if(igs->second == gsfIs[iEle]) {
00482 isGoodBrem = true;
00483 sortedBrem_index = temp_sortedBrem_index;
00484 }
00485 enteredInPrimaryGsf = true;
00486 break;
00487 }
00488 }
00489 if(enteredInPrimaryGsf)
00490 break;
00491 }
00492
00493 if(isGoodBrem) {
00494
00495
00496
00497 std::multimap<double, unsigned int> ecalOtherGsfElems;
00498 block.associatedElements( ecalBrem_index,linkData,
00499 ecalOtherGsfElems,
00500 reco::PFBlockElement::GSF,
00501 reco::PFBlock::LINKTEST_ALL);
00502 if (ecalOtherGsfElems.size() > 0) {
00503 const reco::PFBlockElementGsfTrack * gsfCheck =
00504 dynamic_cast<const reco::PFBlockElementGsfTrack*>((&elements[ecalOtherGsfElems.begin()->second]));
00505 if(ecalOtherGsfElems.begin()->second != gsfIs[iEle] &&
00506 gsfCheck->trackType(reco::PFBlockElement::T_FROM_GAMMACONV) == false) {
00507 continue;
00508 }
00509 }
00510
00511 const reco::PFBlockElementBrem * BremEl =
00512 dynamic_cast<const reco::PFBlockElementBrem*>((&elements[sortedBrem_index]));
00513
00514 reco::PFClusterRef clusterRef =
00515 elements[ecalBrem_index].clusterRef();
00516
00517
00518 float sortedBremEcal_deta = fabs(clusterRef->position().eta() - BremEl->positionAtECALEntrance().eta());
00519
00520
00521 if(sortedBremEcal_deta < 0.015) {
00522
00523 cleanedEcalBremElems.insert(pair<unsigned int,unsigned int>(sortedBrem_index,ecalBrem_index));
00524
00525 unsigned int BremTrajP = BremEl->indTrajPoint();
00526 if (BremTrajP > latestBrem_trajP) {
00527 latestBrem_trajP = BremTrajP;
00528 latestBrem_index = sortedBrem_index;
00529 }
00530 if (DebugSetLinksDetailed)
00531 cout << " brem Index " << sortedBrem_index
00532 << " associated cluster " << ecalBrem_index << " BremTrajP " << BremTrajP <<endl;
00533
00534
00535
00536
00537 localactive[ecalBrem_index] = false;
00538 bool alreadyfound = false;
00539 for(unsigned int ii=0;ii<keyBremIndex.size();ii++) {
00540 if (sortedBrem_index == keyBremIndex[ii]) alreadyfound = true;
00541 }
00542 if (alreadyfound == false) {
00543 keyBremIndex.push_back(sortedBrem_index);
00544 localactive[sortedBrem_index] = false;
00545 }
00546 }
00547 }
00548 }
00549 }
00550
00551
00552
00553 vector<unsigned int> GsfElemIndex(0);
00554 vector<unsigned int> EcalIndex(0);
00555
00556
00557 if (ClosestEcalGsf_index < CutIndex) {
00558 GsfElemIndex.push_back(ClosestEcalGsf_index);
00559 localactive[ClosestEcalGsf_index] = false;
00560 for (std::multimap<double, unsigned int>::iterator ii = ecalGsfElems.begin();
00561 ii != ecalGsfElems.end();ii++) {
00562 if(localactive[ii->second]) {
00563
00564 std::multimap<double, unsigned int> ecalOtherGsfElems;
00565 block.associatedElements( ii->second,linkData,
00566 ecalOtherGsfElems,
00567 reco::PFBlockElement::GSF,
00568 reco::PFBlock::LINKTEST_ALL);
00569 if(ecalOtherGsfElems.size()) {
00570 if(ecalOtherGsfElems.begin()->second != gsfIs[iEle]) continue;
00571 }
00572
00573
00574 reco::PFClusterRef clusterRef = elements[(ii->second)].clusterRef();
00575 float etacl = clusterRef->eta();
00576 if( fabs(eta_gsf-etacl) < 0.05) {
00577 GsfElemIndex.push_back(ii->second);
00578 localactive[ii->second] = false;
00579 if (DebugSetLinksDetailed)
00580 cout << " ExtraCluster From Gsf " << ii->second << endl;
00581 }
00582 }
00583 }
00584 }
00585
00586
00587
00588
00589
00590
00591
00592
00593
00594
00595
00596
00597
00598
00599
00600
00601
00602
00603
00604
00605
00606
00607
00608
00609
00610
00611
00612
00613
00614
00615
00616
00617 if(GsfElemIndex.size() == 0){
00618 if(latestBrem_index < CutIndex) {
00619 unsigned int ckey = cleanedEcalBremElems.count(latestBrem_index);
00620 if(ckey == 1) {
00621 unsigned int temp_cal =
00622 cleanedEcalBremElems.find(latestBrem_index)->second;
00623 GsfElemIndex.push_back(temp_cal);
00624 if (DebugSetLinksDetailed)
00625 cout << "******************** Gsf Cluster From Brem " << temp_cal
00626 << " Latest Brem index " << latestBrem_index
00627 << " ************************* " << endl;
00628 }
00629 else{
00630 pair<multimap<unsigned int,unsigned int>::iterator,multimap<unsigned int,unsigned int>::iterator> ret;
00631 ret = cleanedEcalBremElems.equal_range(latestBrem_index);
00632 multimap<unsigned int,unsigned int>::iterator it;
00633 for(it=ret.first; it!=ret.second; ++it) {
00634 GsfElemIndex.push_back((*it).second);
00635 if (DebugSetLinksDetailed)
00636 cout << "******************** Gsf Cluster From Brem " << (*it).second
00637 << " Latest Brem index " << latestBrem_index
00638 << " ************************* " << endl;
00639 }
00640 }
00641
00642 unsigned int elToErase = 0;
00643 for(unsigned int i = 0; i<keyBremIndex.size();i++) {
00644 if(latestBrem_index == keyBremIndex[i]) {
00645 elToErase = i;
00646 }
00647 }
00648 keyBremIndex.erase(keyBremIndex.begin()+elToErase);
00649 }
00650 }
00651
00652
00653
00654
00655
00656 for(unsigned int iConv=0; iConv<gsfIs.size(); iConv++) {
00657 if(iConv != iEle) {
00658
00659 const reco::PFBlockElementGsfTrack * gsfConv =
00660 dynamic_cast<const reco::PFBlockElementGsfTrack*>((&elements[(gsfIs[iConv])]));
00661
00662
00663 if(gsfConv->trackType(reco::PFBlockElement::T_FROM_GAMMACONV)){
00664 if (DebugSetLinksDetailed)
00665 cout << " PFElectronAlgo:: I'm running on convGsfBrem " << endl;
00666
00667 float conv_dist = block.dist(gsfIs[iConv],gsfIs[iEle],
00668 linkData,reco::PFBlock::LINKTEST_ALL);
00669 if(conv_dist > 0.) {
00670
00671
00672 std::multimap<double, unsigned int> ecalConvElems;
00673 block.associatedElements( gsfIs[iConv],linkData,
00674 ecalConvElems,
00675 reco::PFBlockElement::ECAL,
00676 reco::PFBlock::LINKTEST_ALL );
00677 if(ecalConvElems.size() > 0) {
00678
00679 if(localactive[(ecalConvElems.begin()->second)] == true) {
00680 if (DebugSetLinksDetailed)
00681 cout << " PFElectronAlgo:: convGsfBrem has a ECAL cluster linked and free" << endl;
00682
00683 std::multimap<double, unsigned int> ecalOtherGsfPrimElems;
00684 block.associatedElements( ecalConvElems.begin()->second,linkData,
00685 ecalOtherGsfPrimElems,
00686 reco::PFBlockElement::GSF,
00687 reco::PFBlock::LINKTEST_ALL);
00688 if(ecalOtherGsfPrimElems.size()>0) {
00689 unsigned int gsfprimcheck_index = ecalOtherGsfPrimElems.begin()->second;
00690 const reco::PFBlockElementGsfTrack * gsfCheck =
00691 dynamic_cast<const reco::PFBlockElementGsfTrack*>((&elements[gsfprimcheck_index]));
00692 if(gsfCheck->trackType(reco::PFBlockElement::T_FROM_GAMMACONV) == false) continue;
00693
00694 reco::PFClusterRef clusterRef = elements[ecalConvElems.begin()->second].clusterRef();
00695 if (DebugSetLinksDetailed)
00696 cout << " PFElectronAlgo: !!!!!!! convGsfBrem ECAL cluster has been stored !!!!!!! "
00697 << " Energy " << clusterRef->energy() << " eta,phi " << clusterRef->position().eta()
00698 <<", " << clusterRef->position().phi() << endl;
00699
00700 GsfElemIndex.push_back(ecalConvElems.begin()->second);
00701 convGsfTrack_.push_back(make_pair(ecalConvElems.begin()->second,gsfIs[iConv]));
00702 localactive[ecalConvElems.begin()->second] = false;
00703
00704 }
00705 }
00706 }
00707 }
00708 }
00709 }
00710 }
00711
00712
00713
00714 EcalIndex.insert(EcalIndex.end(),GsfElemIndex.begin(),GsfElemIndex.end());
00715
00716
00717
00718
00719 for(unsigned int i =0;i<keyBremIndex.size();i++) {
00720 unsigned int ikey = keyBremIndex[i];
00721 unsigned int ckey = cleanedEcalBremElems.count(ikey);
00722 vector<unsigned int> BremElemIndex(0);
00723 if(ckey == 1) {
00724 unsigned int temp_cal =
00725 cleanedEcalBremElems.find(ikey)->second;
00726 BremElemIndex.push_back(temp_cal);
00727 }
00728 else{
00729 pair<multimap<unsigned int,unsigned int>::iterator,multimap<unsigned int,unsigned int>::iterator> ret;
00730 ret = cleanedEcalBremElems.equal_range(ikey);
00731 multimap<unsigned int,unsigned int>::iterator it;
00732 for(it=ret.first; it!=ret.second; ++it) {
00733 BremElemIndex.push_back((*it).second);
00734 }
00735 }
00736 EcalIndex.insert(EcalIndex.end(),BremElemIndex.begin(),BremElemIndex.end());
00737 associatedToBrems_.insert(pair<unsigned int,vector<unsigned int> >(ikey,BremElemIndex));
00738 }
00739
00740
00741
00742 vector<unsigned int> convBremKFTrack;
00743 convBremKFTrack.clear();
00744 if (kfElems.size() > 0) {
00745 for(std::multimap<double, unsigned int>::iterator itkf = kfElems.begin();
00746 itkf != kfElems.end(); ++itkf) {
00747 const reco::PFBlockElementTrack * TrkEl =
00748 dynamic_cast<const reco::PFBlockElementTrack*>((&elements[itkf->second]));
00749 bool isPrim = isPrimaryTrack(*TrkEl,*GsfEl);
00750
00751 if(!isPrim) {
00752
00753
00754 std::multimap<double, unsigned int> ecalConvElems;
00755 block.associatedElements( itkf->second,linkData,
00756 ecalConvElems,
00757 reco::PFBlockElement::ECAL,
00758 reco::PFBlock::LINKTEST_ALL );
00759 if(ecalConvElems.size() > 0) {
00760
00761 TrackRef trkRef = TrkEl->trackRef();
00762
00763 unsigned int Algo = whichTrackAlgo(trkRef);
00764
00765 float secpin = trkRef->p();
00766
00767 const reco::PFBlockElementCluster * clust =
00768 dynamic_cast<const reco::PFBlockElementCluster*>((&elements[(ecalConvElems.begin()->second)]));
00769 float eneclust =clust->clusterRef()->energy();
00770
00771
00772
00773
00774
00775 std::multimap<double, unsigned int> hcalConvElems;
00776 block.associatedElements( itkf->second,linkData,
00777 hcalConvElems,
00778 reco::PFBlockElement::HCAL,
00779 reco::PFBlock::LINKTEST_ALL );
00780
00781 bool isHoHE = false;
00782 bool isHoE = false;
00783 bool isPoHE = false;
00784
00785 float enehcalclust = -1;
00786 if(hcalConvElems.size() > 0) {
00787 const reco::PFBlockElementCluster * clusthcal =
00788 dynamic_cast<const reco::PFBlockElementCluster*>((&elements[(hcalConvElems.begin()->second)]));
00789 enehcalclust =clusthcal->clusterRef()->energy();
00790
00791 if( (enehcalclust / (enehcalclust+eneclust) ) > 0.1 && Algo < 3) {
00792 isHoHE = true;
00793 if(enehcalclust > eneclust)
00794 isHoE = true;
00795 if(secpin > (enehcalclust+eneclust) )
00796 isPoHE = true;
00797 }
00798 }
00799
00800
00801 if(localactive[(ecalConvElems.begin()->second)] == false) {
00802
00803 if(isHoE || isPoHE) {
00804 if (DebugSetLinksDetailed)
00805 cout << "PFElectronAlgo:: LOCKED ECAL REJECTED TRACK FOR H/E or P/(H+E) "
00806 << " H/H+E " << enehcalclust/(enehcalclust+eneclust)
00807 << " H/E " << enehcalclust/eneclust
00808 << " P/(H+E) " << secpin/(enehcalclust+eneclust)
00809 << " HCAL ENE " << enehcalclust
00810 << " ECAL ENE " << eneclust
00811 << " secPIN " << secpin
00812 << " Algo Track " << Algo << endl;
00813 continue;
00814 }
00815
00816
00817 for(unsigned int iecal =0; iecal < EcalIndex.size(); iecal++) {
00818
00819
00820 if(EcalIndex[iecal] == ecalConvElems.begin()->second) {
00821 if (DebugSetLinksDetailed)
00822 cout << " PFElectronAlgo:: Conv Brem Recovery locked cluster and I will lock also the KF track " << endl;
00823 convBremKFTrack.push_back(itkf->second);
00824 }
00825 }
00826 }
00827 else{
00828
00829
00830
00831 if(isHoHE){
00832 if (DebugSetLinksDetailed)
00833 cout << "PFElectronAlgo:: FREE ECAL REJECTED TRACK FOR H/H+E "
00834 << " H/H+E " << (enehcalclust / (enehcalclust+eneclust) )
00835 << " H/E " << enehcalclust/eneclust
00836 << " P/(H+E) " << secpin/(enehcalclust+eneclust)
00837 << " HCAL ENE " << enehcalclust
00838 << " ECAL ENE " << eneclust
00839 << " secPIN " << secpin
00840 << " Algo Track " << Algo << endl;
00841 continue;
00842 }
00843
00844
00845 std::multimap<double, unsigned int> ecalOtherKFPrimElems;
00846 block.associatedElements( ecalConvElems.begin()->second,linkData,
00847 ecalOtherKFPrimElems,
00848 reco::PFBlockElement::TRACK,
00849 reco::PFBlock::LINKTEST_ALL);
00850 if(ecalOtherKFPrimElems.size() > 0) {
00851
00852
00853
00854 bool isFromGSF = false;
00855 for(std::multimap<double, unsigned int>::iterator itclos = kfElems.begin();
00856 itclos != kfElems.end(); ++itclos) {
00857 if(ecalOtherKFPrimElems.begin()->second == itclos->second) {
00858 isFromGSF = true;
00859 break;
00860 }
00861 }
00862 if(isFromGSF){
00863
00864
00865
00866
00867 float Epin = eneclust/secpin;
00868
00869
00870 float totenergy = 0.;
00871 for(unsigned int ikeyecal = 0;
00872 ikeyecal<EcalIndex.size(); ikeyecal++){
00873
00874 bool foundcluster = false;
00875 if(ikeyecal > 0) {
00876 for(unsigned int i2 = 0; i2<ikeyecal-1; i2++) {
00877 if(EcalIndex[ikeyecal] == EcalIndex[i2])
00878 foundcluster = true;
00879 }
00880 }
00881 if(foundcluster) continue;
00882 const reco::PFBlockElementCluster * clusasso =
00883 dynamic_cast<const reco::PFBlockElementCluster*>((&elements[(EcalIndex[ikeyecal])]));
00884 totenergy += clusasso->clusterRef()->energy();
00885 }
00886
00887
00888
00889 if(totenergy == 0.) {
00890 if (DebugSetLinksDetailed)
00891 cout << "PFElectronAlgo:: REJECTED_NULLTOT totenergy " << totenergy << endl;
00892 continue;
00893 }
00894
00895
00896 if(Epin > 3) {
00897 double res_before = fabs((totenergy-Pin_gsf)/Pin_gsf);
00898 double res_after = fabs(((totenergy+eneclust)-Pin_gsf)/Pin_gsf);
00899
00900 if(res_before < res_after) {
00901 if (DebugSetLinksDetailed)
00902 cout << "PFElectronAlgo::REJECTED_RES totenergy " << totenergy << " Pin_gsf " << Pin_gsf << " cluster to secondary " << eneclust
00903 << " Res before " << res_before << " res_after " << res_after << endl;
00904 continue;
00905 }
00906 }
00907
00908 if (DebugSetLinksDetailed)
00909 cout << "PFElectronAlgo:: conv brem found asso to ECAL linked to a secondary KF " << endl;
00910 convBremKFTrack.push_back(itkf->second);
00911 GsfElemIndex.push_back(ecalConvElems.begin()->second);
00912 EcalIndex.push_back(ecalConvElems.begin()->second);
00913 localactive[(ecalConvElems.begin()->second)] = false;
00914 localactive[(itkf->second)] = false;
00915 }
00916 }
00917 }
00918 }
00919 }
00920 }
00921 }
00922
00923
00924 if(EcalIndex.size() > 0 && useEGammaSupercluster_) {
00925 double sumEtEcalInTheCone = 0.;
00926
00927
00928 const reco::PFBlockElementCluster * clust =
00929 dynamic_cast<const reco::PFBlockElementCluster*>((&elements[EcalIndex[0]]));
00930 double PhiFC = clust->clusterRef()->position().Phi();
00931 double EtaFC = clust->clusterRef()->position().Eta();
00932
00933
00934 for(unsigned int iEcal=0; iEcal<ecalIs.size(); iEcal++) {
00935 bool foundcluster = false;
00936 for(unsigned int ikeyecal = 0;
00937 ikeyecal<EcalIndex.size(); ikeyecal++){
00938 if(ecalIs[iEcal] == EcalIndex[ikeyecal])
00939 foundcluster = true;
00940 }
00941
00942
00943 if(foundcluster == false) {
00944 const reco::PFBlockElementCluster * clustExt =
00945 dynamic_cast<const reco::PFBlockElementCluster*>((&elements[ecalIs[iEcal]]));
00946 double eta_clust = clustExt->clusterRef()->position().Eta();
00947 double phi_clust = clustExt->clusterRef()->position().Phi();
00948 double theta_clust = clustExt->clusterRef()->position().Theta();
00949 double deta_clust = eta_clust - EtaFC;
00950 double dphi_clust = phi_clust - PhiFC;
00951 if ( dphi_clust < -M_PI )
00952 dphi_clust = dphi_clust + 2.*M_PI;
00953 else if ( dphi_clust > M_PI )
00954 dphi_clust = dphi_clust - 2.*M_PI;
00955 double DR = sqrt(deta_clust*deta_clust+
00956 dphi_clust*dphi_clust);
00957
00958
00959 if(fabs(deta_clust) > 0.05 && DR < coneEcalIsoForEgammaSC_) {
00960 vector<double> ps1Ene(0);
00961 vector<double> ps2Ene(0);
00962 double ps1,ps2;
00963 ps1=ps2=0.;
00964 double EE_calib = pfcalib_.energyEm(*(clustExt->clusterRef()),ps1Ene,ps2Ene,ps1,ps2,applyCrackCorrections_);
00965 double ET_calib = EE_calib*sin(theta_clust);
00966 sumEtEcalInTheCone += ET_calib;
00967 }
00968 }
00969 }
00970
00971
00972 unsigned int sumNTracksInTheCone = 0;
00973 double sumPtTracksInTheCone = 0.;
00974 for(unsigned int iTrack=0; iTrack<trackIs.size(); iTrack++) {
00975
00976 if(localactive[(trackIs[iTrack])]==true) {
00977 const reco::PFBlockElementTrack * kfEle =
00978 dynamic_cast<const reco::PFBlockElementTrack*>((&elements[(trackIs[iTrack])]));
00979 reco::TrackRef trkref = kfEle->trackRef();
00980 if (trkref.isNonnull()) {
00981 double deta_trk = trkref->eta() - RefGSF->etaMode();
00982 double dphi_trk = trkref->phi() - RefGSF->phiMode();
00983 if ( dphi_trk < -M_PI )
00984 dphi_trk = dphi_trk + 2.*M_PI;
00985 else if ( dphi_trk > M_PI )
00986 dphi_trk = dphi_trk - 2.*M_PI;
00987 double DR = sqrt(deta_trk*deta_trk+
00988 dphi_trk*dphi_trk);
00989
00990 reco::HitPattern kfHitPattern = trkref->hitPattern();
00991 int NValPixelHit = kfHitPattern.numberOfValidPixelHits();
00992
00993 if(DR < coneTrackIsoForEgammaSC_ && NValPixelHit >=3) {
00994 sumNTracksInTheCone++;
00995 sumPtTracksInTheCone+=trkref->pt();
00996 }
00997 }
00998 }
00999 }
01000
01001
01002 bool isBarrelIsolated = false;
01003 if( fabs(EtaFC < 1.478) &&
01004 (sumEtEcalInTheCone < sumEtEcalIsoForEgammaSC_barrel_ &&
01005 (sumNTracksInTheCone < nTrackIsoForEgammaSC_ || sumPtTracksInTheCone < sumPtTrackIsoForEgammaSC_barrel_)))
01006 isBarrelIsolated = true;
01007
01008
01009 bool isEndcapIsolated = false;
01010 if( fabs(EtaFC >= 1.478) &&
01011 (sumEtEcalInTheCone < sumEtEcalIsoForEgammaSC_endcap_ &&
01012 (sumNTracksInTheCone < nTrackIsoForEgammaSC_ || sumPtTracksInTheCone < sumPtTrackIsoForEgammaSC_endcap_)))
01013 isEndcapIsolated = true;
01014
01015
01016
01017 if(DebugSetLinksDetailed) {
01018 if(fabs(EtaFC < 1.478) && isBarrelIsolated == false) {
01019 cout << "**** PFElectronAlgo:: SUPERCLUSTER FOUND BUT FAILS ISOLATION:BARREL *** "
01020 << " sumEtEcalInTheCone " <<sumEtEcalInTheCone
01021 << " sumNTracksInTheCone " << sumNTracksInTheCone
01022 << " sumPtTracksInTheCone " << sumPtTracksInTheCone << endl;
01023 }
01024 if(fabs(EtaFC >= 1.478) && isEndcapIsolated == false) {
01025 cout << "**** PFElectronAlgo:: SUPERCLUSTER FOUND BUT FAILS ISOLATION:ENDCAP *** "
01026 << " sumEtEcalInTheCone " <<sumEtEcalInTheCone
01027 << " sumNTracksInTheCone " << sumNTracksInTheCone
01028 << " sumPtTracksInTheCone " << sumPtTracksInTheCone << endl;
01029 }
01030 }
01031
01032
01033
01034
01035 if(isBarrelIsolated || isEndcapIsolated ) {
01036
01037
01038
01039 double totenergy = 0.;
01040 for(unsigned int ikeyecal = 0;
01041 ikeyecal<EcalIndex.size(); ikeyecal++){
01042
01043 bool foundcluster = false;
01044 if(ikeyecal > 0) {
01045 for(unsigned int i2 = 0; i2<ikeyecal-1; i2++) {
01046 if(EcalIndex[ikeyecal] == EcalIndex[i2])
01047 foundcluster = true;;
01048 }
01049 }
01050 if(foundcluster) continue;
01051 const reco::PFBlockElementCluster * clusasso =
01052 dynamic_cast<const reco::PFBlockElementCluster*>((&elements[(EcalIndex[ikeyecal])]));
01053 totenergy += clusasso->clusterRef()->energy();
01054 }
01055
01056
01057
01058
01059 for(unsigned int ikeyecal = 0;
01060 ikeyecal<EcalIndex.size(); ikeyecal++){
01061
01062 bool foundcluster = false;
01063 if(ikeyecal > 0) {
01064 for(unsigned int i2 = 0; i2<ikeyecal-1; i2++) {
01065 if(EcalIndex[ikeyecal] == EcalIndex[i2])
01066 foundcluster = true;
01067 }
01068 }
01069 if(foundcluster) continue;
01070
01071
01072 std::multimap<double, unsigned int> ecalFromSuperClusterElems;
01073 block.associatedElements( EcalIndex[ikeyecal],linkData,
01074 ecalFromSuperClusterElems,
01075 reco::PFBlockElement::ECAL,
01076 reco::PFBlock::LINKTEST_ALL);
01077 if(ecalFromSuperClusterElems.size() > 0) {
01078 for(std::multimap<double, unsigned int>::iterator itsc = ecalFromSuperClusterElems.begin();
01079 itsc != ecalFromSuperClusterElems.end(); ++itsc) {
01080 if(localactive[itsc->second] == false) {
01081 continue;
01082 }
01083
01084 std::multimap<double, unsigned int> ecalOtherKFPrimElems;
01085 block.associatedElements( itsc->second,linkData,
01086 ecalOtherKFPrimElems,
01087 reco::PFBlockElement::TRACK,
01088 reco::PFBlock::LINKTEST_ALL);
01089 if(ecalOtherKFPrimElems.size() > 0) {
01090 if(localactive[ecalOtherKFPrimElems.begin()->second] == true) {
01091 if (DebugSetLinksDetailed)
01092 cout << "**** PFElectronAlgo:: SUPERCLUSTER FOUND BUT FAILS KF VETO *** " << endl;
01093 continue;
01094 }
01095 }
01096 bool isInTheEtaRange = false;
01097 const reco::PFBlockElementCluster * clustToAdd =
01098 dynamic_cast<const reco::PFBlockElementCluster*>((&elements[itsc->second]));
01099 double deta_clustToAdd = clustToAdd->clusterRef()->position().Eta() - EtaFC;
01100 double ene_clustToAdd = clustToAdd->clusterRef()->energy();
01101
01102 if(fabs(deta_clustToAdd) < 0.05)
01103 isInTheEtaRange = true;
01104
01105
01106 bool isBetterEpin = false;
01107 if(isInTheEtaRange == false ) {
01108 if (DebugSetLinksDetailed)
01109 cout << "**** PFElectronAlgo:: SUPERCLUSTER FOUND BUT FAILS GAMMA DETA RANGE *** "
01110 << fabs(deta_clustToAdd) << endl;
01111
01112 if(KfGsf_index < CutIndex) {
01113
01114 double res_before_gsf = fabs((totenergy-Pin_gsf)/Pin_gsf);
01115 double res_after_gsf = fabs(((totenergy+ene_clustToAdd)-Pin_gsf)/Pin_gsf);
01116
01117
01118 const reco::PFBlockElementTrack * trackEl =
01119 dynamic_cast<const reco::PFBlockElementTrack*>((&elements[KfGsf_index]));
01120 double Pin_kf = trackEl->trackRef()->p();
01121 double res_before_kf = fabs((totenergy-Pin_kf)/Pin_kf);
01122 double res_after_kf = fabs(((totenergy+ene_clustToAdd)-Pin_kf)/Pin_kf);
01123
01124
01125 if(res_after_gsf < res_before_gsf && res_after_kf < res_before_kf ) {
01126 isBetterEpin = true;
01127 }
01128 else {
01129 if (DebugSetLinksDetailed)
01130 cout << "**** PFElectronAlgo:: SUPERCLUSTER FOUND AND FAILS ALSO RES_EPIN"
01131 << " tot energy " << totenergy
01132 << " Pin_gsf " << Pin_gsf
01133 << " Pin_kf " << Pin_kf
01134 << " cluster from SC to ADD " << ene_clustToAdd
01135 << " Res before GSF " << res_before_gsf << " res_after_gsf " << res_after_gsf
01136 << " Res before KF " << res_before_kf << " res_after_kf " << res_after_kf << endl;
01137 }
01138 }
01139 }
01140
01141 if(isInTheEtaRange || isBetterEpin) {
01142 if (DebugSetLinksDetailed)
01143 cout << "!!!! PFElectronAlgo:: ECAL from SUPERCLUSTER FOUND !!!!! " << endl;
01144 GsfElemIndex.push_back(itsc->second);
01145 EcalIndex.push_back(itsc->second);
01146 localactive[(itsc->second)] = false;
01147 }
01148 }
01149 }
01150 }
01151 }
01152 }
01153
01154
01155 if(KfGsf_index < CutIndex)
01156 GsfElemIndex.push_back(KfGsf_index);
01157 else if(KfGsf_secondIndex < CutIndex)
01158 GsfElemIndex.push_back(KfGsf_secondIndex);
01159
01160
01161 GsfElemIndex.insert(GsfElemIndex.end(),convBremKFTrack.begin(),convBremKFTrack.end());
01162 GsfElemIndex.insert(GsfElemIndex.end(),keyBremIndex.begin(),keyBremIndex.end());
01163 associatedToGsf_.insert(pair<unsigned int, vector<unsigned int> >(gsfIs[iEle],GsfElemIndex));
01164
01165
01166 for(unsigned int ikeyecal = 0;
01167 ikeyecal<EcalIndex.size(); ikeyecal++){
01168
01169
01170 vector<unsigned int> EcalElemsIndex(0);
01171
01172 std::multimap<double, unsigned int> PS1Elems;
01173 block.associatedElements( EcalIndex[ikeyecal],linkData,
01174 PS1Elems,
01175 reco::PFBlockElement::PS1,
01176 reco::PFBlock::LINKTEST_ALL );
01177 for( std::multimap<double, unsigned int>::iterator it = PS1Elems.begin();
01178 it != PS1Elems.end();it++) {
01179 unsigned int index = it->second;
01180 if(localactive[index] == true) {
01181
01182
01183 std::multimap<double, unsigned> sortedECAL;
01184 block.associatedElements( index, linkData,
01185 sortedECAL,
01186 reco::PFBlockElement::ECAL,
01187 reco::PFBlock::LINKTEST_ALL );
01188 unsigned jEcal = sortedECAL.begin()->second;
01189 if ( jEcal != EcalIndex[ikeyecal]) continue;
01190
01191
01192 EcalElemsIndex.push_back(index);
01193 localactive[index] = false;
01194 }
01195 }
01196
01197 std::multimap<double, unsigned int> PS2Elems;
01198 block.associatedElements( EcalIndex[ikeyecal],linkData,
01199 PS2Elems,
01200 reco::PFBlockElement::PS2,
01201 reco::PFBlock::LINKTEST_ALL );
01202 for( std::multimap<double, unsigned int>::iterator it = PS2Elems.begin();
01203 it != PS2Elems.end();it++) {
01204 unsigned int index = it->second;
01205 if(localactive[index] == true) {
01206
01207 std::multimap<double, unsigned> sortedECAL;
01208 block.associatedElements( index, linkData,
01209 sortedECAL,
01210 reco::PFBlockElement::ECAL,
01211 reco::PFBlock::LINKTEST_ALL );
01212 unsigned jEcal = sortedECAL.begin()->second;
01213 if ( jEcal != EcalIndex[ikeyecal]) continue;
01214
01215 EcalElemsIndex.push_back(index);
01216 localactive[index] = false;
01217 }
01218 }
01219 if(ikeyecal == 0) {
01220
01221
01222
01223 std::multimap<double, unsigned int> hcalGsfElems;
01224 block.associatedElements( gsfIs[iEle],linkData,
01225 hcalGsfElems,
01226 reco::PFBlockElement::HCAL,
01227 reco::PFBlock::LINKTEST_ALL );
01228 for( std::multimap<double, unsigned int>::iterator it = hcalGsfElems.begin();
01229 it != hcalGsfElems.end();it++) {
01230 unsigned int index = it->second;
01231
01232
01233
01234
01235
01236
01237
01238
01239
01240
01241
01242
01243 EcalElemsIndex.push_back(index);
01244 localactive[index] = false;
01245
01246
01247 }
01248
01249 if(hcalGsfElems.size() == 0 && ClosestEcalWithKf == true) {
01250 std::multimap<double, unsigned int> hcalKfElems;
01251 block.associatedElements( KfGsf_index,linkData,
01252 hcalKfElems,
01253 reco::PFBlockElement::HCAL,
01254 reco::PFBlock::LINKTEST_ALL );
01255 for( std::multimap<double, unsigned int>::iterator it = hcalKfElems.begin();
01256 it != hcalKfElems.end();it++) {
01257 unsigned int index = it->second;
01258 if(localactive[index] == true) {
01259
01260
01261 std::multimap<double, unsigned> sortedKf;
01262 block.associatedElements( index, linkData,
01263 sortedKf,
01264 reco::PFBlockElement::TRACK,
01265 reco::PFBlock::LINKTEST_ALL );
01266 unsigned jKf = sortedKf.begin()->second;
01267 if ( jKf != KfGsf_index) continue;
01268 EcalElemsIndex.push_back(index);
01269 localactive[index] = false;
01270 }
01271 }
01272 }
01273
01274 std::multimap<double, unsigned int> kfEtraElems;
01275 block.associatedElements( EcalIndex[ikeyecal],linkData,
01276 kfEtraElems,
01277 reco::PFBlockElement::TRACK,
01278 reco::PFBlock::LINKTEST_ALL );
01279 if(kfEtraElems.size() > 0) {
01280 for( std::multimap<double, unsigned int>::iterator it = kfEtraElems.begin();
01281 it != kfEtraElems.end();it++) {
01282 unsigned int index = it->second;
01283
01284
01285
01286
01287
01288
01289
01290 bool thisIsAMuon = false;
01291 thisIsAMuon = PFMuonAlgo::isMuon(elements[index]);
01292 if (DebugSetLinksDetailed && thisIsAMuon)
01293 cout << " This is a Muon: index " << index << endl;
01294 if(localactive[index] == true && !thisIsAMuon) {
01295 if(index != KfGsf_index) {
01296
01297
01298 std::multimap<double, unsigned> sortedECAL;
01299 block.associatedElements( index, linkData,
01300 sortedECAL,
01301 reco::PFBlockElement::ECAL,
01302 reco::PFBlock::LINKTEST_ALL );
01303 unsigned jEcal = sortedECAL.begin()->second;
01304 if ( jEcal != EcalIndex[ikeyecal]) continue;
01305 EcalElemsIndex.push_back(index);
01306 localactive[index] = false;
01307 }
01308 }
01309 }
01310 }
01311
01312 }
01313 associatedToEcal_.insert(pair<unsigned int,vector<unsigned int> >(EcalIndex[ikeyecal],EcalElemsIndex));
01314 }
01315 }
01316 }
01317
01318
01319
01320
01321
01322 if (DebugSetLinksSummary) {
01323 if(IsThereAGoodGSFTrack) {
01324 if (DebugSetLinksSummary) cout << " -- The Link Summary --" << endl;
01325 for(map<unsigned int,vector<unsigned int> >::iterator it = associatedToGsf_.begin();
01326 it != associatedToGsf_.end(); it++) {
01327
01328 if (DebugSetLinksSummary) cout << " AssoGsf " << it->first << endl;
01329 vector<unsigned int> eleasso = it->second;
01330 for(unsigned int i=0;i<eleasso.size();i++) {
01331 PFBlockElement::Type type = elements[eleasso[i]].type();
01332 if(type == reco::PFBlockElement::BREM) {
01333 if (DebugSetLinksSummary)
01334 cout << " AssoGsfElements BREM " << eleasso[i] << endl;
01335 }
01336 else if(type == reco::PFBlockElement::ECAL) {
01337 if (DebugSetLinksSummary)
01338 cout << " AssoGsfElements ECAL " << eleasso[i] << endl;
01339 }
01340 else if(type == reco::PFBlockElement::TRACK) {
01341 if (DebugSetLinksSummary)
01342 cout << " AssoGsfElements KF " << eleasso[i] << endl;
01343 }
01344 else {
01345 if (DebugSetLinksSummary)
01346 cout << " AssoGsfElements ????? " << eleasso[i] << endl;
01347 }
01348 }
01349 }
01350
01351 for(map<unsigned int,vector<unsigned int> >::iterator it = associatedToBrems_.begin();
01352 it != associatedToBrems_.end(); it++) {
01353 if (DebugSetLinksSummary) cout << " AssoBrem " << it->first << endl;
01354 vector<unsigned int> eleasso = it->second;
01355 for(unsigned int i=0;i<eleasso.size();i++) {
01356 PFBlockElement::Type type = elements[eleasso[i]].type();
01357 if(type == reco::PFBlockElement::ECAL) {
01358 if (DebugSetLinksSummary)
01359 cout << " AssoBremElements ECAL " << eleasso[i] << endl;
01360 }
01361 else {
01362 if (DebugSetLinksSummary)
01363 cout << " AssoBremElements ????? " << eleasso[i] << endl;
01364 }
01365 }
01366 }
01367
01368 for(map<unsigned int,vector<unsigned int> >::iterator it = associatedToEcal_.begin();
01369 it != associatedToEcal_.end(); it++) {
01370 if (DebugSetLinksSummary) cout << " AssoECAL " << it->first << endl;
01371 vector<unsigned int> eleasso = it->second;
01372 for(unsigned int i=0;i<eleasso.size();i++) {
01373 PFBlockElement::Type type = elements[eleasso[i]].type();
01374 if(type == reco::PFBlockElement::PS1) {
01375 if (DebugSetLinksSummary)
01376 cout << " AssoECALElements PS1 " << eleasso[i] << endl;
01377 }
01378 else if(type == reco::PFBlockElement::PS2) {
01379 if (DebugSetLinksSummary)
01380 cout << " AssoECALElements PS2 " << eleasso[i] << endl;
01381 }
01382 else if(type == reco::PFBlockElement::HCAL) {
01383 if (DebugSetLinksSummary)
01384 cout << " AssoECALElements HCAL " << eleasso[i] << endl;
01385 }
01386 else {
01387 if (DebugSetLinksSummary)
01388 cout << " AssoHCALElements ????? " << eleasso[i] << endl;
01389 }
01390 }
01391 }
01392 if (DebugSetLinksSummary)
01393 cout << "-- End Summary --" << endl;
01394 }
01395
01396 }
01397
01398 return IsThereAGoodGSFTrack;
01399 }
01400
01401 void PFElectronAlgo::SetIDOutputs(const reco::PFBlockRef& blockRef,
01402 AssMap& associatedToGsf_,
01403 AssMap& associatedToBrems_,
01404 AssMap& associatedToEcal_){
01405 PFEnergyCalibration pfcalib_;
01406 const reco::PFBlock& block = *blockRef;
01407 PFBlock::LinkData linkData = block.linkData();
01408 const edm::OwnVector< reco::PFBlockElement >& elements = block.elements();
01409 bool DebugIDOutputs = false;
01410 if(DebugIDOutputs) cout << " ######## Enter in SetIDOutputs #########" << endl;
01411
01412 unsigned int cgsf=0;
01413 for (map<unsigned int,vector<unsigned int> >::iterator igsf = associatedToGsf_.begin();
01414 igsf != associatedToGsf_.end(); igsf++,cgsf++) {
01415
01416 float Ene_ecalgsf = 0.;
01417 float Ene_hcalgsf = 0.;
01418 double sigmaEtaEta = 0.;
01419 float deta_gsfecal = 0.;
01420 float Ene_ecalbrem = 0.;
01421 float Ene_extraecalgsf = 0.;
01422 bool LateBrem = false;
01423 bool EarlyBrem = false;
01424 int FirstBrem = 1000;
01425 unsigned int ecalGsf_index = 100000;
01426 unsigned int kf_index = 100000;
01427 unsigned int nhits_gsf = 0;
01428 int NumBrem = 0;
01429 reco::TrackRef RefKF;
01430 double posX=0.;
01431 double posY=0.;
01432 double posZ=0.;
01433
01434 unsigned int gsf_index = igsf->first;
01435 const reco::PFBlockElementGsfTrack * GsfEl =
01436 dynamic_cast<const reco::PFBlockElementGsfTrack*>((&elements[gsf_index]));
01437 reco::GsfTrackRef RefGSF = GsfEl->GsftrackRef();
01438 float Ein_gsf = 0.;
01439 if (RefGSF.isNonnull()) {
01440 float m_el=0.00051;
01441 Ein_gsf =sqrt(RefGSF->pMode()*
01442 RefGSF->pMode()+m_el*m_el);
01443 nhits_gsf = RefGSF->hitPattern().trackerLayersWithMeasurement();
01444 }
01445 float Eout_gsf = GsfEl->Pout().t();
01446 float Etaout_gsf = GsfEl->positionAtECALEntrance().eta();
01447
01448
01449 if (DebugIDOutputs)
01450 cout << " setIdOutput! GSF Track: Ein " << Ein_gsf
01451 << " eta,phi " << Etaout_gsf
01452 <<", " << GsfEl->positionAtECALEntrance().phi() << endl;
01453
01454
01455 vector<unsigned int> assogsf_index = igsf->second;
01456 bool FirstEcalGsf = true;
01457 for (unsigned int ielegsf=0;ielegsf<assogsf_index.size();ielegsf++) {
01458 PFBlockElement::Type assoele_type = elements[(assogsf_index[ielegsf])].type();
01459
01460
01461
01462 if(assoele_type == reco::PFBlockElement::TRACK) {
01463 const reco::PFBlockElementTrack * KfTk =
01464 dynamic_cast<const reco::PFBlockElementTrack*>((&elements[(assogsf_index[ielegsf])]));
01465
01466
01467 bool isPrim = isPrimaryTrack(*KfTk,*GsfEl);
01468 if(!isPrim) continue;
01469 RefKF = KfTk->trackRef();
01470 kf_index = assogsf_index[ielegsf];
01471 }
01472
01473
01474 if (assoele_type == reco::PFBlockElement::ECAL) {
01475 unsigned int keyecalgsf = assogsf_index[ielegsf];
01476 vector<unsigned int> assoecalgsf_index = associatedToEcal_.find(keyecalgsf)->second;
01477
01478 vector<double> ps1Ene(0);
01479 vector<double> ps2Ene(0);
01480 for(unsigned int ips =0; ips<assoecalgsf_index.size();ips++) {
01481 PFBlockElement::Type typeassoecal = elements[(assoecalgsf_index[ips])].type();
01482 if (typeassoecal == reco::PFBlockElement::PS1) {
01483 PFClusterRef psref = elements[(assoecalgsf_index[ips])].clusterRef();
01484 ps1Ene.push_back(psref->energy());
01485 }
01486 if (typeassoecal == reco::PFBlockElement::PS2) {
01487 PFClusterRef psref = elements[(assoecalgsf_index[ips])].clusterRef();
01488 ps2Ene.push_back(psref->energy());
01489 }
01490 if (typeassoecal == reco::PFBlockElement::HCAL) {
01491 const reco::PFBlockElementCluster * clust =
01492 dynamic_cast<const reco::PFBlockElementCluster*>((&elements[(assoecalgsf_index[ips])]));
01493 Ene_hcalgsf+=clust->clusterRef()->energy();
01494 }
01495 }
01496 if (FirstEcalGsf) {
01497 FirstEcalGsf = false;
01498 ecalGsf_index = assogsf_index[ielegsf];
01499 reco::PFClusterRef clusterRef = elements[(assogsf_index[ielegsf])].clusterRef();
01500 double ps1,ps2;
01501 ps1=ps2=0.;
01502
01503 Ene_ecalgsf = pfcalib_.energyEm(*clusterRef,ps1Ene,ps2Ene,ps1,ps2,applyCrackCorrections_);
01504
01505 posX += Ene_ecalgsf * clusterRef->position().X();
01506 posY += Ene_ecalgsf * clusterRef->position().Y();
01507 posZ += Ene_ecalgsf * clusterRef->position().Z();
01508
01509 if (DebugIDOutputs)
01510 cout << " setIdOutput! GSF ECAL Cluster E " << Ene_ecalgsf
01511 << " eta,phi " << clusterRef->position().eta()
01512 <<", " << clusterRef->position().phi() << endl;
01513
01514 deta_gsfecal = clusterRef->position().eta() - Etaout_gsf;
01515
01516 vector< const reco::PFCluster * > pfClust_vec(0);
01517 pfClust_vec.clear();
01518 pfClust_vec.push_back(&(*clusterRef));
01519
01520 PFClusterWidthAlgo pfwidth(pfClust_vec);
01521 sigmaEtaEta = pfwidth.pflowSigmaEtaEta();
01522
01523
01524 }
01525 else {
01526 reco::PFClusterRef clusterRef = elements[(assogsf_index[ielegsf])].clusterRef();
01527 float TempClus_energy = pfcalib_.energyEm(*clusterRef,ps1Ene,ps2Ene,applyCrackCorrections_);
01528 Ene_extraecalgsf += TempClus_energy;
01529 posX += TempClus_energy * clusterRef->position().X();
01530 posY += TempClus_energy * clusterRef->position().Y();
01531 posZ += TempClus_energy * clusterRef->position().Z();
01532
01533 if (DebugIDOutputs)
01534 cout << " setIdOutput! Extra ECAL Cluster E "
01535 << TempClus_energy << " Tot " << Ene_extraecalgsf << endl;
01536 }
01537 }
01538
01539
01540 if (assoele_type == reco::PFBlockElement::BREM) {
01541 unsigned int brem_index = assogsf_index[ielegsf];
01542 const reco::PFBlockElementBrem * BremEl =
01543 dynamic_cast<const reco::PFBlockElementBrem*>((&elements[brem_index]));
01544 int TrajPos = (BremEl->indTrajPoint())-2;
01545 if (TrajPos <= 3) EarlyBrem = true;
01546 if (TrajPos < FirstBrem) FirstBrem = TrajPos;
01547
01548 vector<unsigned int> assobrem_index = associatedToBrems_.find(brem_index)->second;
01549 for (unsigned int ibrem = 0; ibrem < assobrem_index.size(); ibrem++){
01550 if (elements[(assobrem_index[ibrem])].type() == reco::PFBlockElement::ECAL) {
01551 unsigned int keyecalbrem = assobrem_index[ibrem];
01552 vector<unsigned int> assoelebrem_index = associatedToEcal_.find(keyecalbrem)->second;
01553 vector<double> ps1EneFromBrem(0);
01554 vector<double> ps2EneFromBrem(0);
01555 for (unsigned int ielebrem=0; ielebrem<assoelebrem_index.size();ielebrem++) {
01556 if (elements[(assoelebrem_index[ielebrem])].type() == reco::PFBlockElement::PS1) {
01557 PFClusterRef psref = elements[(assoelebrem_index[ielebrem])].clusterRef();
01558 ps1EneFromBrem.push_back(psref->energy());
01559 }
01560 if (elements[(assoelebrem_index[ielebrem])].type() == reco::PFBlockElement::PS2) {
01561 PFClusterRef psref = elements[(assoelebrem_index[ielebrem])].clusterRef();
01562 ps2EneFromBrem.push_back(psref->energy());
01563 }
01564 }
01565
01566 if( assobrem_index[ibrem] != ecalGsf_index) {
01567 reco::PFClusterRef clusterRef =
01568 elements[(assobrem_index[ibrem])].clusterRef();
01569 float BremClus_energy = pfcalib_.energyEm(*clusterRef,ps1EneFromBrem,ps2EneFromBrem,applyCrackCorrections_);
01570 Ene_ecalbrem += BremClus_energy;
01571 posX += BremClus_energy * clusterRef->position().X();
01572 posY += BremClus_energy * clusterRef->position().Y();
01573 posZ += BremClus_energy * clusterRef->position().Z();
01574
01575 NumBrem++;
01576 if (DebugIDOutputs) cout << " setIdOutput::BREM Cluster "
01577 << BremClus_energy << " eta,phi "
01578 << clusterRef->position().eta() <<", "
01579 << clusterRef->position().phi() << endl;
01580 }
01581 else {
01582 LateBrem = true;
01583 }
01584 }
01585 }
01586 }
01587 }
01588 if (Ene_ecalgsf > 0.) {
01589
01590
01591
01592 if(RefGSF.isNonnull()) {
01593 PFCandidateElectronExtra myExtra(RefGSF) ;
01594 myExtra.setGsfTrackPout(GsfEl->Pout());
01595 myExtra.setKfTrackRef(RefKF);
01596 float Pt_gsf = RefGSF->ptMode();
01597 lnPt_gsf = log(Pt_gsf);
01598 Eta_gsf = RefGSF->etaMode();
01599
01600
01601 if(RefGSF->ptModeError() > 0.)
01602 dPtOverPt_gsf = RefGSF->ptModeError()/Pt_gsf;
01603
01604 nhit_gsf= RefGSF->hitPattern().trackerLayersWithMeasurement();
01605 chi2_gsf = RefGSF->normalizedChi2();
01606
01607 DPtOverPt_gsf = (RefGSF->ptMode() - GsfEl->Pout().pt())/RefGSF->ptMode();
01608
01609
01610 nhit_kf = 0;
01611 chi2_kf = -0.01;
01612 DPtOverPt_kf = -0.01;
01613 if (RefKF.isNonnull()) {
01614 nhit_kf= RefKF->hitPattern().trackerLayersWithMeasurement();
01615 chi2_kf = RefKF->normalizedChi2();
01616
01617
01618
01619
01620
01621 }
01622
01623
01624 EtotPinMode = (Ene_ecalgsf + Ene_ecalbrem + Ene_extraecalgsf) / Ein_gsf;
01625 EGsfPoutMode = Ene_ecalgsf/Eout_gsf;
01626 EtotBremPinPoutMode = Ene_ecalbrem /(Ein_gsf - Eout_gsf);
01627 DEtaGsfEcalClust = fabs(deta_gsfecal);
01628 myExtra.setSigmaEtaEta(sigmaEtaEta);
01629 myExtra.setDeltaEta(DEtaGsfEcalClust);
01630 SigmaEtaEta = log(sigmaEtaEta);
01631
01632 lateBrem = -1;
01633 firstBrem = -1;
01634 earlyBrem = -1;
01635 if(NumBrem > 0) {
01636 if (LateBrem) lateBrem = 1;
01637 else lateBrem = 0;
01638 firstBrem = FirstBrem;
01639 if(FirstBrem < 4) earlyBrem = 1;
01640 else earlyBrem = 0;
01641 }
01642
01643 HOverHE = Ene_hcalgsf/(Ene_hcalgsf + Ene_ecalgsf);
01644 HOverPin = Ene_hcalgsf / Ein_gsf;
01645 myExtra.setHadEnergy(Ene_hcalgsf);
01646 myExtra.setEarlyBrem(earlyBrem);
01647 myExtra.setLateBrem(lateBrem);
01648
01649
01650
01651 if(DPtOverPt_gsf < -0.2) DPtOverPt_gsf = -0.2;
01652 if(DPtOverPt_gsf > 1.) DPtOverPt_gsf = 1.;
01653
01654 if(dPtOverPt_gsf > 0.3) dPtOverPt_gsf = 0.3;
01655
01656 if(chi2_gsf > 10.) chi2_gsf = 10.;
01657
01658 if(DPtOverPt_kf < -0.2) DPtOverPt_kf = -0.2;
01659 if(DPtOverPt_kf > 1.) DPtOverPt_kf = 1.;
01660
01661 if(chi2_kf > 10.) chi2_kf = 10.;
01662
01663 if(EtotPinMode < 0.) EtotPinMode = 0.;
01664 if(EtotPinMode > 5.) EtotPinMode = 5.;
01665
01666 if(EGsfPoutMode < 0.) EGsfPoutMode = 0.;
01667 if(EGsfPoutMode > 5.) EGsfPoutMode = 5.;
01668
01669 if(EtotBremPinPoutMode < 0.) EtotBremPinPoutMode = 0.01;
01670 if(EtotBremPinPoutMode > 5.) EtotBremPinPoutMode = 5.;
01671
01672 if(DEtaGsfEcalClust > 0.1) DEtaGsfEcalClust = 0.1;
01673
01674 if(SigmaEtaEta < -14) SigmaEtaEta = -14;
01675
01676 if(HOverPin < 0.) HOverPin = 0.;
01677 if(HOverPin > 5.) HOverPin = 5.;
01678 double mvaValue = tmvaReader_->EvaluateMVA("BDT");
01679
01680
01681 BDToutput_[cgsf] = mvaValue;
01682 myExtra.setMVA(mvaValue);
01683 electronExtra_.push_back(myExtra);
01684
01685
01686 if(mvaValue > mvaEleCut_) {
01687
01688
01689
01690
01691
01692
01693
01694
01695
01696
01697
01698
01699
01700 unsigned int iextratrack = 0;
01701 unsigned int itrackHcalLinked = 0;
01702 float SumExtraKfP = 0.;
01703
01704
01705 double Etotal = Ene_ecalgsf + Ene_ecalbrem + Ene_extraecalgsf;
01706 posX /=Etotal;
01707 posY /=Etotal;
01708 posZ /=Etotal;
01709 math::XYZPoint sc_pflow(posX,posY,posZ);
01710 double ETtotal = Etotal*sin(sc_pflow.Theta());
01711 double phiTrack = RefGSF->phiMode();
01712 double dphi_normalsc = sc_pflow.Phi() - phiTrack;
01713 if ( dphi_normalsc < -M_PI )
01714 dphi_normalsc = dphi_normalsc + 2.*M_PI;
01715 else if ( dphi_normalsc > M_PI )
01716 dphi_normalsc = dphi_normalsc - 2.*M_PI;
01717 dphi_normalsc = fabs(dphi_normalsc);
01718
01719
01720 if(ecalGsf_index < 100000) {
01721 vector<unsigned int> assoecalgsf_index = associatedToEcal_.find(ecalGsf_index)->second;
01722 for(unsigned int itrk =0; itrk<assoecalgsf_index.size();itrk++) {
01723 PFBlockElement::Type typeassoecal = elements[(assoecalgsf_index[itrk])].type();
01724 if(typeassoecal == reco::PFBlockElement::TRACK) {
01725 const reco::PFBlockElementTrack * kfTk =
01726 dynamic_cast<const reco::PFBlockElementTrack*>((&elements[(assoecalgsf_index[itrk])]));
01727
01728
01729
01730
01731
01732
01733 reco::TrackRef trackref = kfTk->trackRef();
01734 unsigned int Algo = whichTrackAlgo(trackref);
01735
01736
01737 if(Algo < 3) {
01738 if(DebugIDOutputs)
01739 cout << " The ecalGsf cluster is not isolated: >0 KF extra with algo < 3" << endl;
01740
01741 float p_trk = trackref->p();
01742 SumExtraKfP += p_trk;
01743 iextratrack++;
01744
01745 std::multimap<double, unsigned int> hcalKfElems;
01746 block.associatedElements( assoecalgsf_index[itrk],linkData,
01747 hcalKfElems,
01748 reco::PFBlockElement::HCAL,
01749 reco::PFBlock::LINKTEST_ALL );
01750 if(hcalKfElems.size() > 0) {
01751 itrackHcalLinked++;
01752 }
01753 }
01754 }
01755 }
01756 }
01757 if( iextratrack > 0) {
01758 if(iextratrack > 3 || HOverHE > 0.05 || (SumExtraKfP/Ene_ecalgsf) > 1.
01759 || (ETtotal > 50. && iextratrack > 1 && (Ene_hcalgsf/Ene_ecalgsf) > 0.1) ) {
01760
01761 if(DebugIDOutputs)
01762 cout << " *****This electron candidate is discarded: Non isolated # tracks "
01763 << iextratrack << " HOverHE " << HOverHE
01764 << " SumExtraKfP/Ene_ecalgsf " << SumExtraKfP/Ene_ecalgsf
01765 << " ETtotal " << ETtotal
01766 << " Ene_hcalgsf/Ene_ecalgsf " << Ene_hcalgsf/Ene_ecalgsf
01767 << endl;
01768
01769 BDToutput_[cgsf] = mvaValue-2.;
01770 lockExtraKf_[cgsf] = false;
01771 }
01772
01773
01774 if( (fabs(1.-EtotPinMode) < 0.2 && (fabs(Eta_gsf) < 1.0 || fabs(Eta_gsf) > 2.0)) ||
01775 ((EtotPinMode < 1.1 && EtotPinMode > 0.6) && (fabs(Eta_gsf) >= 1.0 && fabs(Eta_gsf) <= 2.0))) {
01776 if( fabs(1.-EGsfPoutMode) < 0.5 &&
01777 (itrackHcalLinked == iextratrack) &&
01778 kf_index < 100000 ) {
01779
01780 BDToutput_[cgsf] = mvaValue;
01781 lockExtraKf_[cgsf] = false;
01782 if(DebugIDOutputs)
01783 cout << " *****This electron is reactivated # tracks "
01784 << iextratrack << " #tracks hcal linked " << itrackHcalLinked
01785 << " SumExtraKfP/Ene_ecalgsf " << SumExtraKfP/Ene_ecalgsf
01786 << " EtotPinMode " << EtotPinMode << " EGsfPoutMode " << EGsfPoutMode
01787 << " eta gsf " << fabs(Eta_gsf) << " kf index " << kf_index <<endl;
01788 }
01789 }
01790 }
01791
01792 if (HOverPin > 1. && HOverHE > 0.1 && EtotPinMode < 0.5) {
01793 if(DebugIDOutputs)
01794 cout << " *****This electron candidate is discarded HCAL ENERGY "
01795 << " HOverPin " << HOverPin << " HOverHE " << HOverHE << " EtotPinMode" << EtotPinMode << endl;
01796 BDToutput_[cgsf] = mvaValue-4.;
01797 lockExtraKf_[cgsf] = false;
01798 }
01799
01800
01801
01802
01803 if( EtotPinMode < 0.2 && EGsfPoutMode < 0.2 ) {
01804 if(DebugIDOutputs)
01805 cout << " *****This electron candidate is discarded Low ETOTPIN "
01806 << " EtotPinMode " << EtotPinMode << " EGsfPoutMode " << EGsfPoutMode << endl;
01807 BDToutput_[cgsf] = mvaValue-6.;
01808 }
01809
01810
01811 if(ETtotal > 50. && dphi_normalsc > 0.1 ) {
01812 if(DebugIDOutputs)
01813 cout << " *****This electron candidate is discarded Large ANGLE "
01814 << " ETtotal " << ETtotal << " EGsfPoutMode " << dphi_normalsc << endl;
01815 BDToutput_[cgsf] = mvaValue-6.;
01816 }
01817 }
01818
01819
01820 if (DebugIDOutputs) {
01821 cout << " **** BDT observables ****" << endl;
01822 cout << " < Normalization > " << endl;
01823 cout << " lnPt_gsf " << lnPt_gsf
01824 << " Eta_gsf " << Eta_gsf << endl;
01825 cout << " < PureTracking > " << endl;
01826 cout << " dPtOverPt_gsf " << dPtOverPt_gsf
01827 << " DPtOverPt_gsf " << DPtOverPt_gsf
01828 << " chi2_gsf " << chi2_gsf
01829 << " nhit_gsf " << nhit_gsf
01830 << " DPtOverPt_kf " << DPtOverPt_kf
01831 << " chi2_kf " << chi2_kf
01832 << " nhit_kf " << nhit_kf << endl;
01833 cout << " < track-ecal-hcal-ps " << endl;
01834 cout << " EtotPinMode " << EtotPinMode
01835 << " EGsfPoutMode " << EGsfPoutMode
01836 << " EtotBremPinPoutMode " << EtotBremPinPoutMode
01837 << " DEtaGsfEcalClust " << DEtaGsfEcalClust
01838 << " SigmaEtaEta " << SigmaEtaEta
01839 << " HOverHE " << HOverHE
01840 << " HOverPin " << HOverPin
01841 << " lateBrem " << lateBrem
01842 << " firstBrem " << firstBrem << endl;
01843 cout << " !!!!!!!!!!!!!!!! the BDT output !!!!!!!!!!!!!!!!!: direct " << mvaValue
01844 << " corrected " << BDToutput_[cgsf] << endl;
01845
01846 }
01847 }
01848 else {
01849 if (DebugIDOutputs)
01850 cout << " Gsf Ref isNULL " << endl;
01851 BDToutput_[cgsf] = -2.;
01852 }
01853 } else {
01854 if (DebugIDOutputs)
01855 cout << " No clusters associated to the gsf " << endl;
01856 BDToutput_[cgsf] = -2.;
01857 }
01858 }
01859 return;
01860 }
01861
01862
01863
01864 void PFElectronAlgo::SetCandidates(const reco::PFBlockRef& blockRef,
01865 AssMap& associatedToGsf_,
01866 AssMap& associatedToBrems_,
01867 AssMap& associatedToEcal_){
01868
01869 const reco::PFBlock& block = *blockRef;
01870 PFBlock::LinkData linkData = block.linkData();
01871 const edm::OwnVector< reco::PFBlockElement >& elements = block.elements();
01872 PFEnergyResolution pfresol_;
01873 PFEnergyCalibration pfcalib_;
01874
01875 bool DebugIDCandidates = false;
01876
01877
01878
01879 unsigned int cgsf=0;
01880 for (map<unsigned int,vector<unsigned int> >::iterator igsf = associatedToGsf_.begin();
01881 igsf != associatedToGsf_.end(); igsf++,cgsf++) {
01882 unsigned int gsf_index = igsf->first;
01883
01884
01885
01886
01887 int eecal=0;
01888 int hcal=0;
01889 int charge =0;
01890 bool goodphi=true;
01891 math::XYZTLorentzVector momentum_kf,momentum_gsf,momentum,momentum_mean;
01892 float dpt=0; float dpt_kf=0; float dpt_gsf=0; float dpt_mean=0;
01893 float Eene=0; float dene=0; float Hene=0.;
01894 float RawEene = 0.;
01895 double posX=0.;
01896 double posY=0.;
01897 double posZ=0.;
01898 std::vector<float> bremEnergyVec;
01899
01900 float de_gs = 0., de_me = 0., de_kf = 0.;
01901 float m_el=0.00051;
01902 int nhit_kf=0; int nhit_gsf=0;
01903 bool has_gsf=false;
01904 bool has_kf=false;
01905 math::XYZTLorentzVector newmomentum;
01906 float ps1TotEne = 0;
01907 float ps2TotEne = 0;
01908 vector<unsigned int> elementsToAdd(0);
01909 reco::TrackRef RefKF;
01910
01911
01912
01913 elementsToAdd.push_back(gsf_index);
01914 const reco::PFBlockElementGsfTrack * GsfEl =
01915 dynamic_cast<const reco::PFBlockElementGsfTrack*>((&elements[gsf_index]));
01916 const math::XYZPointF& posGsfEcalEntrance = GsfEl->positionAtECALEntrance();
01917 reco::GsfTrackRef RefGSF = GsfEl->GsftrackRef();
01918 if (RefGSF.isNonnull()) {
01919
01920 has_gsf=true;
01921
01922 charge= RefGSF->chargeMode();
01923 nhit_gsf= RefGSF->hitPattern().trackerLayersWithMeasurement();
01924
01925 momentum_gsf.SetPx(RefGSF->pxMode());
01926 momentum_gsf.SetPy(RefGSF->pyMode());
01927 momentum_gsf.SetPz(RefGSF->pzMode());
01928 float ENE=sqrt(RefGSF->pMode()*
01929 RefGSF->pMode()+m_el*m_el);
01930
01931 if( DebugIDCandidates )
01932 cout << "SetCandidates:: GsfTrackRef: Ene " << ENE
01933 << " charge " << charge << " nhits " << nhit_gsf <<endl;
01934
01935 momentum_gsf.SetE(ENE);
01936 dpt_gsf=RefGSF->ptModeError()*
01937 (RefGSF->pMode()/RefGSF->ptMode());
01938
01939 momentum_mean.SetPx(RefGSF->px());
01940 momentum_mean.SetPy(RefGSF->py());
01941 momentum_mean.SetPz(RefGSF->pz());
01942 float ENEm=sqrt(RefGSF->p()*
01943 RefGSF->p()+m_el*m_el);
01944 momentum_mean.SetE(ENEm);
01945 dpt_mean=RefGSF->ptError()*
01946 (RefGSF->p()/RefGSF->pt());
01947 }
01948 else {
01949 if( DebugIDCandidates )
01950 cout << "SetCandidates:: !!!! NULL GSF Track Ref " << endl;
01951 }
01952
01953
01954 vector<unsigned int> assogsf_index = igsf->second;
01955 unsigned int ecalGsf_index = 100000;
01956 bool FirstEcalGsf = true;
01957 for (unsigned int ielegsf=0;ielegsf<assogsf_index.size();ielegsf++) {
01958 PFBlockElement::Type assoele_type = elements[(assogsf_index[ielegsf])].type();
01959 if (assoele_type == reco::PFBlockElement::TRACK) {
01960 elementsToAdd.push_back((assogsf_index[ielegsf]));
01961 const reco::PFBlockElementTrack * KfTk =
01962 dynamic_cast<const reco::PFBlockElementTrack*>((&elements[(assogsf_index[ielegsf])]));
01963
01964 bool isPrim = isPrimaryTrack(*KfTk,*GsfEl);
01965 if(!isPrim) continue;
01966
01967 RefKF = KfTk->trackRef();
01968 if (RefKF.isNonnull()) {
01969 has_kf = true;
01970 dpt_kf=(RefKF->ptError()*RefKF->ptError());
01971 nhit_kf=RefKF->hitPattern().trackerLayersWithMeasurement();
01972 momentum_kf.SetPx(RefKF->px());
01973 momentum_kf.SetPy(RefKF->py());
01974 momentum_kf.SetPz(RefKF->pz());
01975 float ENE=sqrt(RefKF->p()*RefKF->p()+m_el*m_el);
01976 if( DebugIDCandidates )
01977 cout << "SetCandidates:: KFTrackRef: Ene " << ENE << " nhits " << nhit_kf << endl;
01978
01979 momentum_kf.SetE(ENE);
01980 }
01981 else {
01982 if( DebugIDCandidates )
01983 cout << "SetCandidates:: !!!! NULL KF Track Ref " << endl;
01984 }
01985 }
01986
01987 if (assoele_type == reco::PFBlockElement::ECAL) {
01988 unsigned int keyecalgsf = assogsf_index[ielegsf];
01989 vector<unsigned int> assoecalgsf_index = associatedToEcal_.find(keyecalgsf)->second;
01990 vector<double> ps1Ene(0);
01991 vector<double> ps2Ene(0);
01992
01993
01994
01995 for(unsigned int ips =0; ips<assoecalgsf_index.size();ips++) {
01996 PFBlockElement::Type typeassoecal = elements[(assoecalgsf_index[ips])].type();
01997 if (typeassoecal == reco::PFBlockElement::PS1) {
01998 PFClusterRef psref = elements[(assoecalgsf_index[ips])].clusterRef();
01999 ps1Ene.push_back(psref->energy());
02000 elementsToAdd.push_back((assoecalgsf_index[ips]));
02001 }
02002 if (typeassoecal == reco::PFBlockElement::PS2) {
02003 PFClusterRef psref = elements[(assoecalgsf_index[ips])].clusterRef();
02004 ps2Ene.push_back(psref->energy());
02005 elementsToAdd.push_back((assoecalgsf_index[ips]));
02006 }
02007 if (typeassoecal == reco::PFBlockElement::HCAL) {
02008 const reco::PFBlockElementCluster * clust =
02009 dynamic_cast<const reco::PFBlockElementCluster*>((&elements[(assoecalgsf_index[ips])]));
02010 elementsToAdd.push_back((assoecalgsf_index[ips]));
02011 Hene+=clust->clusterRef()->energy();
02012 hcal++;
02013 }
02014 }
02015 elementsToAdd.push_back((assogsf_index[ielegsf]));
02016
02017
02018 const reco::PFBlockElementCluster * clust =
02019 dynamic_cast<const reco::PFBlockElementCluster*>((&elements[(assogsf_index[ielegsf])]));
02020
02021 eecal++;
02022
02023 const reco::PFCluster& cl(*clust->clusterRef());
02024
02025
02026
02027 double ps1,ps2;
02028 ps1=ps2=0.;
02029
02030 float EE = pfcalib_.energyEm(cl,ps1Ene,ps2Ene,ps1,ps2,applyCrackCorrections_);
02031
02032
02033 float ceta=cl.position().eta();
02034 float cphi=cl.position().phi();
02035
02036 float mphi=-2.97025;
02037 if (ceta<0) mphi+=0.00638;
02038 for (int ip=1; ip<19; ip++){
02039 float df= cphi - (mphi+(ip*6.283185/18));
02040 if (fabs(df)<0.01) goodphi=false;
02041 }
02042 float dE=pfresol_.getEnergyResolutionEm(EE,cl.position().eta());
02043 if( DebugIDCandidates )
02044 cout << "SetCandidates:: EcalCluster: EneNoCalib " << clust->clusterRef()->energy()
02045 << " eta,phi " << ceta << "," << cphi << " Calib " << EE << " dE " << dE <<endl;
02046
02047 bool elecCluster=false;
02048 if (FirstEcalGsf) {
02049 FirstEcalGsf = false;
02050 elecCluster=true;
02051 ecalGsf_index = assogsf_index[ielegsf];
02052
02053 RawEene += EE;
02054 }
02055
02056
02057 math::XYZTLorentzVector clusterMomentum;
02058 math::XYZPoint direction=cl.position()/cl.position().R();
02059 clusterMomentum.SetPxPyPzE(EE*direction.x(),
02060 EE*direction.y(),
02061 EE*direction.z(),
02062 EE);
02063 reco::PFCandidate cluster_Candidate((elecCluster)?charge:0,
02064 clusterMomentum,
02065 (elecCluster)? reco::PFCandidate::e : reco::PFCandidate::gamma);
02066
02067 cluster_Candidate.setPs1Energy(ps1);
02068 cluster_Candidate.setPs2Energy(ps2);
02069 cluster_Candidate.setEcalEnergy(EE);
02070
02071
02072
02073 cluster_Candidate.setRawEcalEnergy(EE-ps1-ps2);
02074 cluster_Candidate.setPositionAtECALEntrance(math::XYZPointF(cl.position()));
02075 cluster_Candidate.addElementInBlock(blockRef,assogsf_index[ielegsf]);
02076
02077 std::map<unsigned int,std::vector<reco::PFCandidate> >::iterator itcheck=
02078 electronConstituents_.find(cgsf);
02079 if(itcheck==electronConstituents_.end())
02080 {
02081
02082 std::vector<reco::PFCandidate> tmpVec;
02083 tmpVec.push_back(cluster_Candidate);
02084 electronConstituents_.insert(std::pair<unsigned int, std::vector<reco::PFCandidate> >
02085 (cgsf,tmpVec));
02086 }
02087 else
02088 {
02089 itcheck->second.push_back(cluster_Candidate);
02090 }
02091
02092 Eene+=EE;
02093 posX += EE * cl.position().X();
02094 posY += EE * cl.position().Y();
02095 posZ += EE * cl.position().Z();
02096 ps1TotEne+=ps1;
02097 ps2TotEne+=ps2;
02098 dene+=dE*dE;
02099 }
02100
02101
02102
02103
02104 if (assoele_type == reco::PFBlockElement::BREM) {
02105 unsigned int brem_index = assogsf_index[ielegsf];
02106 vector<unsigned int> assobrem_index = associatedToBrems_.find(brem_index)->second;
02107 elementsToAdd.push_back(brem_index);
02108 for (unsigned int ibrem = 0; ibrem < assobrem_index.size(); ibrem++){
02109 if (elements[(assobrem_index[ibrem])].type() == reco::PFBlockElement::ECAL) {
02110
02111 unsigned int keyecalbrem = assobrem_index[ibrem];
02112 const vector<unsigned int>& assoelebrem_index = associatedToEcal_.find(keyecalbrem)->second;
02113 vector<double> ps1EneFromBrem(0);
02114 vector<double> ps2EneFromBrem(0);
02115 float ps1EneFromBremTot=0.;
02116 float ps2EneFromBremTot=0.;
02117 for (unsigned int ielebrem=0; ielebrem<assoelebrem_index.size();ielebrem++) {
02118 if (elements[(assoelebrem_index[ielebrem])].type() == reco::PFBlockElement::PS1) {
02119 PFClusterRef psref = elements[(assoelebrem_index[ielebrem])].clusterRef();
02120 ps1EneFromBrem.push_back(psref->energy());
02121 ps1EneFromBremTot+=psref->energy();
02122 elementsToAdd.push_back(assoelebrem_index[ielebrem]);
02123 }
02124 if (elements[(assoelebrem_index[ielebrem])].type() == reco::PFBlockElement::PS2) {
02125 PFClusterRef psref = elements[(assoelebrem_index[ielebrem])].clusterRef();
02126 ps2EneFromBrem.push_back(psref->energy());
02127 ps2EneFromBremTot+=psref->energy();
02128 elementsToAdd.push_back(assoelebrem_index[ielebrem]);
02129 }
02130 }
02131
02132
02133 if( assobrem_index[ibrem] != ecalGsf_index) {
02134 elementsToAdd.push_back(assobrem_index[ibrem]);
02135 reco::PFClusterRef clusterRef = elements[(assobrem_index[ibrem])].clusterRef();
02136
02137
02138 double ps1=0;
02139 double ps2=0;
02140 float EE = pfcalib_.energyEm(*clusterRef,ps1EneFromBrem,ps2EneFromBrem,ps1,ps2,applyCrackCorrections_);
02141 bremEnergyVec.push_back(EE);
02142
02143 float ceta = clusterRef->position().eta();
02144
02145 float dE=pfresol_.getEnergyResolutionEm(EE,ceta);
02146 if( DebugIDCandidates )
02147 cout << "SetCandidates:: BremCluster: Ene " << EE << " dE " << dE <<endl;
02148
02149 Eene+=EE;
02150 posX += EE * clusterRef->position().X();
02151 posY += EE * clusterRef->position().Y();
02152 posZ += EE * clusterRef->position().Z();
02153 ps1TotEne+=ps1;
02154 ps2TotEne+=ps2;
02155
02156
02157 dene+=dE*dE;
02158
02159
02160
02161 math::XYZTLorentzVector photonMomentum;
02162 math::XYZPoint direction=clusterRef->position()/clusterRef->position().R();
02163
02164 photonMomentum.SetPxPyPzE(EE*direction.x(),
02165 EE*direction.y(),
02166 EE*direction.z(),
02167 EE);
02168 reco::PFCandidate photon_Candidate(0,photonMomentum, reco::PFCandidate::gamma);
02169
02170 photon_Candidate.setPs1Energy(ps1);
02171 photon_Candidate.setPs2Energy(ps2);
02172 photon_Candidate.setEcalEnergy(EE);
02173
02174
02175
02176 photon_Candidate.setRawEcalEnergy(EE-ps1-ps2);
02177 photon_Candidate.setPositionAtECALEntrance(math::XYZPointF(clusterRef->position()));
02178 photon_Candidate.addElementInBlock(blockRef,assobrem_index[ibrem]);
02179
02180
02181 std::map<unsigned int,std::vector<reco::PFCandidate> >::iterator itcheck=
02182 electronConstituents_.find(cgsf);
02183 if(itcheck==electronConstituents_.end())
02184 {
02185
02186 std::vector<reco::PFCandidate> tmpVec;
02187 tmpVec.push_back(photon_Candidate);
02188 electronConstituents_.insert(std::pair<unsigned int, std::vector<reco::PFCandidate> >
02189 (cgsf,tmpVec));
02190 }
02191 else
02192 {
02193 itcheck->second.push_back(photon_Candidate);
02194 }
02195 }
02196 }
02197 }
02198 }
02199 }
02200 if (has_gsf) {
02201
02202
02203 double unCorrEene = Eene;
02204 double absEta = fabs(momentum_gsf.Eta());
02205 double emTheta = momentum_gsf.Theta();
02206 if( DebugIDCandidates )
02207 cout << "PFEelectronAlgo:: absEta " << absEta << " theta " << emTheta
02208 << " EneRaw " << Eene << " Err " << dene;
02209
02210
02211 if(usePFSCEleCalib_ && (unCorrEene*sin(emTheta)) < 200 && unCorrEene > 0.) {
02212 if( absEta < 1.5) {
02213 double Etene = Eene*sin(emTheta);
02214 double emCorrFull_et = thePFSCEnergyCalibration_->SCCorrEtEtaBarrel(Etene, absEta);
02215 Eene = emCorrFull_et/sin(emTheta);
02216 }
02217 else {
02218 double Etene = Eene*sin(emTheta);
02219 double emCorrFull_et = thePFSCEnergyCalibration_->SCCorrEtEtaEndcap(Etene, absEta);
02220 Eene = emCorrFull_et/sin(emTheta);
02221 }
02222 dene = sqrt(dene)*(Eene/unCorrEene);
02223 dene = dene*dene;
02224 }
02225
02226 if( DebugIDCandidates )
02227 cout << " EneCorrected " << Eene << " Err " << dene << endl;
02228
02229
02230
02231
02232
02233 if(has_kf && unCorrEene > 0.) {
02234 posX /=unCorrEene;
02235 posY /=unCorrEene;
02236 posZ /=unCorrEene;
02237 math::XYZPoint sc_pflow(posX,posY,posZ);
02238
02239 std::multimap<double, unsigned int> bremElems;
02240 block.associatedElements( gsf_index,linkData,
02241 bremElems,
02242 reco::PFBlockElement::BREM,
02243 reco::PFBlock::LINKTEST_ALL );
02244
02245 double phiTrack = RefGSF->phiMode();
02246 if(bremElems.size()>0) {
02247 unsigned int brem_index = bremElems.begin()->second;
02248 const reco::PFBlockElementBrem * BremEl =
02249 dynamic_cast<const reco::PFBlockElementBrem*>((&elements[brem_index]));
02250 phiTrack = BremEl->positionAtECALEntrance().phi();
02251 }
02252
02253 double dphi_normalsc = sc_pflow.Phi() - phiTrack;
02254 if ( dphi_normalsc < -M_PI )
02255 dphi_normalsc = dphi_normalsc + 2.*M_PI;
02256 else if ( dphi_normalsc > M_PI )
02257 dphi_normalsc = dphi_normalsc - 2.*M_PI;
02258
02259 int chargeGsf = RefGSF->chargeMode();
02260 int chargeKf = RefKF->charge();
02261
02262 int chargeSC = 0;
02263 if(dphi_normalsc < 0.)
02264 chargeSC = 1;
02265 else
02266 chargeSC = -1;
02267
02268 if(chargeKf == chargeGsf)
02269 charge = chargeGsf;
02270 else if(chargeGsf == chargeSC)
02271 charge = chargeGsf;
02272 else
02273 charge = chargeKf;
02274
02275 if( DebugIDCandidates )
02276 cout << "PFElectronAlgo:: charge determination "
02277 << " charge GSF " << chargeGsf
02278 << " charge KF " << chargeKf
02279 << " charge SC " << chargeSC
02280 << " Final Charge " << charge << endl;
02281
02282 }
02283
02284
02285 if ((nhit_gsf<8) && (has_kf)){
02286
02287
02288
02289 momentum=momentum_kf;
02290 float Fe=Eene;
02291 float scale= Fe/momentum.E();
02292
02293
02294 if (Eene < 0.0001) {
02295 Fe = momentum.E();
02296 scale = 1.;
02297 }
02298
02299
02300 newmomentum.SetPxPyPzE(scale*momentum.Px(),
02301 scale*momentum.Py(),
02302 scale*momentum.Pz(),Fe);
02303 if( DebugIDCandidates )
02304 cout << "SetCandidates:: (nhit_gsf<8) && (has_kf):: pt " << newmomentum.pt() << " Ene " << Fe <<endl;
02305
02306
02307 }
02308 if ((nhit_gsf>7) || (has_kf==false)){
02309 if(Eene > 0.0001) {
02310 de_gs=1-momentum_gsf.E()/Eene;
02311 de_me=1-momentum_mean.E()/Eene;
02312 de_kf=1-momentum_kf.E()/Eene;
02313 }
02314
02315 momentum=momentum_gsf;
02316 dpt=1/(dpt_gsf*dpt_gsf);
02317
02318 if(dene > 0.)
02319 dene= 1./dene;
02320
02321 float Fe = 0.;
02322 if(Eene > 0.0001) {
02323 Fe =((dene*Eene) +(dpt*momentum.E()))/(dene+dpt);
02324 }
02325 else {
02326 Fe=momentum.E();
02327 }
02328
02329 if ((de_gs>0.05)&&(de_kf>0.05)){
02330 Fe=Eene;
02331 }
02332 if ((de_gs<-0.1)&&(de_me<-0.1) &&(de_kf<0.) &&
02333 (momentum.E()/dpt_gsf) > 5. && momentum_gsf.pt() < 30.){
02334 Fe=momentum.E();
02335 }
02336 float scale= Fe/momentum.E();
02337
02338 newmomentum.SetPxPyPzE(scale*momentum.Px(),
02339 scale*momentum.Py(),
02340 scale*momentum.Pz(),Fe);
02341 if( DebugIDCandidates )
02342 cout << "SetCandidates::(nhit_gsf>7) || (has_kf==false) " << newmomentum.pt() << " Ene " << Fe <<endl;
02343
02344
02345 }
02346 if (newmomentum.pt()>0.5){
02347
02348
02349
02350
02351 if( DebugIDCandidates )
02352 cout << "SetCandidates:: I am before doing candidate " <<endl;
02353
02354
02355 std::vector<float> clusterEnergyVec;
02356 clusterEnergyVec.push_back(RawEene);
02357 clusterEnergyVec.insert(clusterEnergyVec.end(),bremEnergyVec.begin(),bremEnergyVec.end());
02358
02359
02360 std::vector<reco::PFCandidateElectronExtra>::iterator itextra;
02361 PFElectronExtraEqual myExtraEqual(RefGSF);
02362 itextra=find_if(electronExtra_.begin(),electronExtra_.end(),myExtraEqual);
02363 if(itextra!=electronExtra_.end()) {
02364 itextra->setClusterEnergies(clusterEnergyVec);
02365 }
02366 else {
02367 if(RawEene>0.)
02368 std::cout << " There is a big problem with the electron extra, PFElectronAlgo should crash soon " << RawEene << std::endl;
02369 }
02370
02371 reco::PFCandidate::ParticleType particleType
02372 = reco::PFCandidate::e;
02373 reco::PFCandidate temp_Candidate;
02374 temp_Candidate = PFCandidate(charge,newmomentum,particleType);
02375 temp_Candidate.set_mva_e_pi(BDToutput_[cgsf]);
02376 temp_Candidate.setEcalEnergy(Eene);
02377 temp_Candidate.setRawEcalEnergy(RawEene);
02378
02379 temp_Candidate.setHcalEnergy(Hene);
02380 temp_Candidate.setPs1Energy(ps1TotEne);
02381 temp_Candidate.setPs2Energy(ps2TotEne);
02382 temp_Candidate.setTrackRef(RefKF);
02383
02384 temp_Candidate.setGsfTrackRef(RefGSF);
02385 temp_Candidate.setPositionAtECALEntrance(posGsfEcalEntrance);
02386
02387 temp_Candidate.setVertex(RefGSF->vertex());
02388
02389
02390 if( DebugIDCandidates )
02391 cout << "SetCandidates:: I am after doing candidate " <<endl;
02392
02393 for (unsigned int elad=0; elad<elementsToAdd.size();elad++){
02394 temp_Candidate.addElementInBlock(blockRef,elementsToAdd[elad]);
02395 }
02396
02397 if(BDToutput_[cgsf] >= mvaEleCut_)
02398 elCandidate_.push_back(temp_Candidate);
02399
02400
02401 reco::PFCandidate extendedElCandidate(temp_Candidate);
02402
02403 std::map<unsigned int, std::vector<reco::PFCandidate> >::const_iterator itcluster=
02404 electronConstituents_.find(cgsf);
02405 if(itcluster!=electronConstituents_.end())
02406 {
02407 const std::vector<reco::PFCandidate> & theClusters=itcluster->second;
02408 unsigned nclus=theClusters.size();
02409
02410 for(unsigned iclus=0;iclus<nclus;++iclus)
02411 {
02412 extendedElCandidate.addDaughter(theClusters[iclus]);
02413 }
02414 }
02415
02416
02417
02418
02419
02420 allElCandidate_.push_back(extendedElCandidate);
02421
02422 }
02423 else {
02424 BDToutput_[cgsf] = -1.;
02425
02426 if( DebugIDCandidates )
02427 cout << "SetCandidates:: No Candidate Produced because of Pt cut: 0.5 " <<endl;
02428 }
02429 }
02430 else {
02431 BDToutput_[cgsf] = -1.;
02432 if( DebugIDCandidates )
02433 cout << "SetCandidates:: No Candidate Produced because of No GSF Track Ref " <<endl;
02434 }
02435 }
02436 return;
02437 }
02438
02439
02440
02441 void PFElectronAlgo::SetActive(const reco::PFBlockRef& blockRef,
02442 AssMap& associatedToGsf_,
02443 AssMap& associatedToBrems_,
02444 AssMap& associatedToEcal_,
02445 std::vector<bool>& active){
02446 const reco::PFBlock& block = *blockRef;
02447 PFBlock::LinkData linkData = block.linkData();
02448
02449 const edm::OwnVector< reco::PFBlockElement >& elements = block.elements();
02450
02451 unsigned int cgsf=0;
02452 for (map<unsigned int,vector<unsigned int> >::iterator igsf = associatedToGsf_.begin();
02453 igsf != associatedToGsf_.end(); igsf++,cgsf++) {
02454
02455
02456 if(BDToutput_[cgsf] < mvaEleCut_) continue;
02457
02458 unsigned int gsf_index = igsf->first;
02459 active[gsf_index] = false;
02460 vector<unsigned int> assogsf_index = igsf->second;
02461 for (unsigned int ielegsf=0;ielegsf<assogsf_index.size();ielegsf++) {
02462 PFBlockElement::Type assoele_type = elements[(assogsf_index[ielegsf])].type();
02463
02464 active[(assogsf_index[ielegsf])] = false;
02465 if (assoele_type == reco::PFBlockElement::ECAL) {
02466 unsigned int keyecalgsf = assogsf_index[ielegsf];
02467
02468
02469 if(fifthStepKfTrack_.size() > 0) {
02470 for(unsigned int itr = 0; itr < fifthStepKfTrack_.size(); itr++) {
02471 if(fifthStepKfTrack_[itr].first == keyecalgsf) {
02472 active[(fifthStepKfTrack_[itr].second)] = false;
02473 }
02474 }
02475 }
02476
02477
02478 if(convGsfTrack_.size() > 0) {
02479 for(unsigned int iconv = 0; iconv < convGsfTrack_.size(); iconv++) {
02480 if(convGsfTrack_[iconv].first == keyecalgsf) {
02481
02482 active[(convGsfTrack_[iconv].second)] = false;
02483
02484 std::multimap<double, unsigned> convKf;
02485 block.associatedElements( convGsfTrack_[iconv].second,
02486 linkData,
02487 convKf,
02488 reco::PFBlockElement::TRACK,
02489 reco::PFBlock::LINKTEST_ALL );
02490 if(convKf.size() > 0) {
02491 active[convKf.begin()->second] = false;
02492 }
02493 }
02494 }
02495 }
02496
02497
02498 vector<unsigned int> assoecalgsf_index = associatedToEcal_.find(keyecalgsf)->second;
02499 for(unsigned int ips =0; ips<assoecalgsf_index.size();ips++) {
02500
02501 if (elements[(assoecalgsf_index[ips])].type() == reco::PFBlockElement::PS1)
02502 active[(assoecalgsf_index[ips])] = false;
02503 if (elements[(assoecalgsf_index[ips])].type() == reco::PFBlockElement::PS2)
02504 active[(assoecalgsf_index[ips])] = false;
02505 if (elements[(assoecalgsf_index[ips])].type() == reco::PFBlockElement::TRACK) {
02506 if(lockExtraKf_[cgsf] == true) {
02507 active[(assoecalgsf_index[ips])] = false;
02508 }
02509 }
02510 }
02511 }
02512 if (assoele_type == reco::PFBlockElement::BREM) {
02513 unsigned int brem_index = assogsf_index[ielegsf];
02514 vector<unsigned int> assobrem_index = associatedToBrems_.find(brem_index)->second;
02515 for (unsigned int ibrem = 0; ibrem < assobrem_index.size(); ibrem++){
02516 if (elements[(assobrem_index[ibrem])].type() == reco::PFBlockElement::ECAL) {
02517 unsigned int keyecalbrem = assobrem_index[ibrem];
02518
02519 active[(assobrem_index[ibrem])] = false;
02520
02521
02522 if(fifthStepKfTrack_.size() > 0) {
02523 for(unsigned int itr = 0; itr < fifthStepKfTrack_.size(); itr++) {
02524 if(fifthStepKfTrack_[itr].first == keyecalbrem) {
02525 active[(fifthStepKfTrack_[itr].second)] = false;
02526 }
02527 }
02528 }
02529
02530 vector<unsigned int> assoelebrem_index = associatedToEcal_.find(keyecalbrem)->second;
02531
02532 for (unsigned int ielebrem=0; ielebrem<assoelebrem_index.size();ielebrem++) {
02533 if (elements[(assoelebrem_index[ielebrem])].type() == reco::PFBlockElement::PS1)
02534 active[(assoelebrem_index[ielebrem])] = false;
02535 if (elements[(assoelebrem_index[ielebrem])].type() == reco::PFBlockElement::PS2)
02536 active[(assoelebrem_index[ielebrem])] = false;
02537 }
02538 }
02539 }
02540 }
02541 }
02542 }
02543 return;
02544 }
02545 void PFElectronAlgo::setEGElectronCollection(const reco::GsfElectronCollection & egelectrons) {
02546 theGsfElectrons_ = & egelectrons;
02547 }
02548 unsigned int PFElectronAlgo::whichTrackAlgo(const reco::TrackRef& trackRef) {
02549 unsigned int Algo = 0;
02550 switch (trackRef->algo()) {
02551 case TrackBase::ctf:
02552 case TrackBase::iter0:
02553 case TrackBase::iter1:
02554 Algo = 0;
02555 break;
02556 case TrackBase::iter2:
02557 Algo = 1;
02558 break;
02559 case TrackBase::iter3:
02560 Algo = 2;
02561 break;
02562 case TrackBase::iter4:
02563 Algo = 3;
02564 break;
02565 case TrackBase::iter5:
02566 Algo = 4;
02567 break;
02568 default:
02569 Algo = 5;
02570 break;
02571 }
02572 return Algo;
02573 }
02574 bool PFElectronAlgo::isPrimaryTrack(const reco::PFBlockElementTrack& KfEl,
02575 const reco::PFBlockElementGsfTrack& GsfEl) {
02576 bool isPrimary = false;
02577
02578 GsfPFRecTrackRef gsfPfRef = GsfEl.GsftrackRefPF();
02579
02580 if(gsfPfRef.isNonnull()) {
02581 PFRecTrackRef kfPfRef = KfEl.trackRefPF();
02582 PFRecTrackRef kfPfRef_fromGsf = (*gsfPfRef).kfPFRecTrackRef();
02583 if(kfPfRef.isNonnull() && kfPfRef_fromGsf.isNonnull()) {
02584 reco::TrackRef kfref= (*kfPfRef).trackRef();
02585 reco::TrackRef kfref_fromGsf = (*kfPfRef_fromGsf).trackRef();
02586 if(kfref.isNonnull() && kfref_fromGsf.isNonnull()) {
02587 if(kfref == kfref_fromGsf)
02588 isPrimary = true;
02589 }
02590 }
02591 }
02592
02593 return isPrimary;
02594 }