CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_4_1_8_patch13/src/RecoParticleFlow/PFProducer/src/PFElectronAlgo.cc

Go to the documentation of this file.
00001 //
00002 // Original Author: Daniele Benedetti: daniele.benedetti@cern.ch
00003 //
00004 
00005 #include "RecoParticleFlow/PFProducer/interface/PFElectronAlgo.h"
00006 #include "RecoParticleFlow/PFProducer/interface/PFMuonAlgo.h"
00007 //#include "DataFormats/ParticleFlowReco/interface/PFBlockElementGsfTrack.h"
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   // Set the tmva reader
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   //tmvaReader_->AddVariable("nhit_gsf",&nhit_gsf);
00064   tmvaReader_->AddVariable("chi2_gsf",&chi2_gsf);
00065   //tmvaReader_->AddVariable("DPtOverPt_kf",&DPtOverPt_kf);
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 //   tmvaReader_->AddVariable("HOverPin",&HOverPin);
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   // the maps are initialized 
00083   AssMap associatedToGsf;
00084   AssMap associatedToBrems;
00085   AssMap associatedToEcal;
00086 
00087   // should be cleaned as often as often as possible
00088   elCandidate_.clear();
00089   electronExtra_.clear();
00090   allElCandidate_.clear();
00091   electronConstituents_.clear();
00092   fifthStepKfTrack_.clear();
00093   convGsfTrack_.clear();
00094   // SetLinks finds all the elements (kf,ecal,ps,hcal,brems) 
00095   // associated to each gsf track
00096   bool blockHasGSF = SetLinks(blockRef,associatedToGsf,
00097                               associatedToBrems,associatedToEcal,
00098                               active);
00099   
00100   // check if there is at least a gsf track in the block. 
00101   if (blockHasGSF) {
00102     
00103     BDToutput_.clear();
00104 
00105     lockExtraKf_.clear();
00106     // For each GSF track is initialized a BDT value = -1
00107     BDToutput_.assign(associatedToGsf.size(),-1.);
00108     lockExtraKf_.assign(associatedToGsf.size(),true);
00109 
00110     // The FinalID is run and BDToutput values is assigned 
00111     SetIDOutputs(blockRef,associatedToGsf,associatedToBrems,associatedToEcal);  
00112 
00113     // For each GSF track that pass the BDT configurable cut a pf candidate electron is created. 
00114     // This function finds also the best estimation of the initial electron 4-momentum.
00115 
00116     SetCandidates(blockRef,associatedToGsf,associatedToBrems,associatedToEcal);
00117     if (elCandidate_.size() > 0 ){
00118       isvalid_ = true;
00119       // when a pfelectron candidate is created all the elements associated to the
00120       // electron are locked. 
00121       SetActive(blockRef,associatedToGsf,associatedToBrems,associatedToEcal,active);
00122     }
00123   } // endif blockHasGSF
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   // no other cut are not used anymore. We use the default of PFBlockAlgo
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   // Save the elements in shorter vectors like in PFAlgo.
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       // Check if the track is already identified as a muon
00160       thisIsAMuon =  PFMuonAlgo::isMuon(elements[iEle]);
00161       // Otherwise store index
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       // Check if the track has a KF partner identified as a muon
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       // Otherwise store index
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   // ******************* Start Link *****************************
00196   // Do something only if a gsf track is found in the block
00197   if(IsThereAGSFTrack) {
00198     
00199 
00200     // LocalLock the Elements associated to a Kf tracks and not to a Gsf
00201     // The clusters associated both to a kf track and to a brem tangend 
00202     // are then assigned only to the kf track
00203     // Could be improved doing this after. 
00204 
00205     // 19 Mar 2010 adding the KF track from Gamma Conv. 
00206     // They are linked to the GSF tracks they are not considered
00207     // anymore in the following ecal cluster locking 
00208 
00209     //     cout<<"#########################################################"<<endl;
00210     //     cout<<"#####           Process Block:                      #####"<<endl;
00211     //     cout<<"#########################################################"<<endl;
00212     //     cout<<block<<endl;
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         // This means that the considered kf is *not* associated
00222         // to any gsf track
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             // Check if this clusters is however well linked to a primary gsf track
00232             // if this the case the cluster is not locked.
00233             
00234             bool isGsfLinked = false;
00235             for(unsigned int iGsf=0; iGsf<gsfIs.size(); iGsf++) {  
00236               // if the ecal cluster is associated contemporary to a KF track
00237               // and to a GSF track from conv, it is assigned to the KF track 
00238               // In this way we can loose some cluster but it is safer for double counting. 
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               // add protection against energy loss because
00256               // of the tracking fifth step
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       } // gsfElems.size()
00273     } // loop on kf tracks
00274 
00275 
00276     // start loop on gsf tracks
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       // if GsfTrack fron converted bremsstralung continue
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       // Find Associated Kf Track elements and Ecal to KF elements
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         // 19 Mar 2010 now a loop is needed because > 1 KF track could
00312         // be associated to the same GSF track
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             // Find clusters associated to kftrack using linkbyrechit
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       // Find the closest Ecal clusters associated to this Gsf
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           // check energy compatibility for outer eta != ecal entrance, looping tracks
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             // Check that this cluster is not closer to another primary Gsf track
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               // get if it is closed to a conv brem gsf tracks
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           // do not lock at the moment we need this for the late brem
00383         }
00384       }
00385       // if any cluster is found with the gsf-ecal link, try with kf-ecal
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           // Check if this cluster is not closer to another Gsf track
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       //  Find the brems associated to this Gsf
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         // Find the ecal clusters associated to the brems
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           //to be changed, using the distance
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           // Find the closest brem
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           // check that this brem is that one coming from the same *primary* gsf
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             //  Check that this cluster is not closer to another Gsf Track
00496             // The check is not performed on KF track because the ecal clusters are aready locked.
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             // Triangular cut on plan chi2:deta -> OLD
00520             //if((0.0075*sortedBremEcal_chi2 + 100.*sortedBremEcal_deta -1.5) < 0.) {
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               // > 1 ecal clusters could be associated to the same brem twice: allowed N-1 link. 
00535               // But the brem need to be stored once. 
00536               // locallock the brem and the ecal clusters
00537               localactive[ecalBrem_index] = false;  // the cluster
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;   // the brem
00545               }
00546             }
00547           }
00548         }
00549       }
00550 
00551       
00552       // Find Possible Extra Cluster associated to the gsf/kf
00553       vector<unsigned int> GsfElemIndex(0);
00554       vector<unsigned int> EcalIndex(0);
00555 
00556       // locallock the ecal cluster associated to the gsf
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             // Check that this cluster is not closer to another Gsf Track
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             // get the cluster only if the deta (ecal-gsf) < 0.05
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       //Add the possibility to link other ecal clusters from kf. 
00587      
00588 //       for (std::multimap<double, unsigned int>::iterator ii = ecalKfElems.begin();
00589 //         ii != ecalKfElems.end();ii++) {
00590 //      if(localactive[ii->second]) {
00591 //         // Check that this cluster is not closer to another Gsf Track    
00592 //        std::multimap<double, unsigned int> ecalOtherGsfElems;
00593 //        block.associatedElements( ii->second,linkData,
00594 //                                  ecalOtherGsfElems,
00595 //                                  reco::PFBlockElement::GSF,
00596 //                                  reco::PFBlock::LINKTEST_CHI2);
00597 //        if(ecalOtherGsfElems.size()) {
00598 //          if(ecalOtherGsfElems.begin()->second != gsfIs[iEle]) continue;
00599 //        } 
00600 //        GsfElemIndex.push_back(ii->second);
00601 //        reco::PFClusterRef clusterRef = elements[(ii->second)].clusterRef();
00602 //        float etacl =  clusterRef->eta();
00603 //        if( fabs(eta_gsf-etacl) < 0.05) {         
00604 //          localactive[ii->second] = false;
00605 //          if (DebugSetLinksDetailed)
00606 //            cout << " ExtraCluster From KF " << ii->second << endl;
00607 //        }
00608 //      }
00609 //       }
00610       
00611       //****************** Fill Maps *************************
00612 
00613       // The GsfMap    
00614 
00615       // if any clusters have been associated to the gsf track    
00616       // use the Ecal clusters associated to the latest brem and associate it to the gsf
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           // erase the brem. 
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       // Get Extra Clusters from converted brem gsf tracks. The locallock method
00653       // tells me if the ecal cluster has been already assigned to the primary
00654       // gsf track or to a brem
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           // look at only to secondary gsf tracks
00663           if(gsfConv->trackType(reco::PFBlockElement::T_FROM_GAMMACONV)){
00664             if (DebugSetLinksDetailed)
00665               cout << "  PFElectronAlgo:: I'm running on convGsfBrem " << endl;
00666             // check if they are linked to the primary
00667             float conv_dist = block.dist(gsfIs[iConv],gsfIs[iEle],
00668                                          linkData,reco::PFBlock::LINKTEST_ALL);
00669             if(conv_dist > 0.) {
00670               // find the closest ecal cluster associated to conversions
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                 // the ecal cluster is still active?
00679                 if(localactive[(ecalConvElems.begin()->second)] == true) {
00680                   if (DebugSetLinksDetailed)
00681                     cout << "  PFElectronAlgo:: convGsfBrem has a ECAL cluster linked and free" << endl;
00682                   // Check that this cluster is not closer to another primary Gsf track
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       // The BremMap
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       // 19 Mar 2010: add KF and ECAL elements from converted brem photons
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             // search for linked ECAL clusters
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               // Further Cleaning: DANIELE This could be improved!
00761               TrackRef trkRef =   TrkEl->trackRef();
00762               // iter0, iter1, iter2, iter3 = Algo < 3
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               //1)  ******* Reject secondary KF tracks linked to also an HCAL cluster with H/(E+H) > 0.1
00772               //            This is applied also to KF linked to locked ECAL cluster
00773               //            NOTE: trusting the H/(E+H) and not the conv brem selection increse the number
00774               //                  of charged hadrons around the electron. DANIELE? re-think about this. 
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                 // NOTE: DANIELE? Are you sure you want to use the Algo type here? 
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                 // check if this track has been alread assigned to an ECAL cluster
00817                 for(unsigned int iecal =0; iecal < EcalIndex.size(); iecal++) {
00818                   // in case this track is already assigned to a locked ECAL cluster
00819                   // the secondary kf track is also saved for further lock
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                 // ECAL cluster free
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                 // check that this cluster is not cluser to another KF track (primary)
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                   // check that this ECAL clusters is the best associated to at least one of the  KF tracks
00853                   // linked to the considered GSF track
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                     // Further Cleaning: DANIELE This could be improved!                                    
00865                   
00866                                           
00867                     float Epin = eneclust/secpin;
00868                     
00869                     // compute the pfsupercluster energy till now
00870                     float totenergy = 0.;
00871                     for(unsigned int ikeyecal = 0; 
00872                         ikeyecal<EcalIndex.size(); ikeyecal++){
00873                       // EcalIndex can have the same cluster save twice (because of the late brem cluster).
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                     // Further Cleaning: DANIELE This could be improved! 
00888                     //2) *****  Do not consider secondary tracks if the GSF and brems have failed in finding ECAL clusters
00889                     if(totenergy == 0.) {
00890                       if (DebugSetLinksDetailed)
00891                         cout << "PFElectronAlgo:: REJECTED_NULLTOT totenergy " << totenergy << endl;
00892                       continue;
00893                     }
00894                     
00895                     //3) ****** Reject secondary KF tracks that have an high E/secPin and that make worse the Etot/pin 
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       // 4May import EG supercluster
00924       if(EcalIndex.size() > 0 && useEGammaSupercluster_) {
00925         double sumEtEcalInTheCone  = 0.;
00926         
00927         // Position of the first cluster
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         // Compute ECAL isolation ->
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           // -> only for clusters not already in the PFSCCluster
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             //Jurassic veto in deta
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         } //EndLoop Additional ECAL clusters in the block
00970         
00971         // Compute track isolation: number of tracks && sumPt
00972         unsigned int sumNTracksInTheCone = 0;
00973         double sumPtTracksInTheCone = 0.;
00974         for(unsigned int iTrack=0; iTrack<trackIs.size(); iTrack++) {
00975           // the track from the electron are already locked at this stage
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         // only print out
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           // Compute TotEnergy
01039           double totenergy = 0.;
01040           for(unsigned int ikeyecal = 0; 
01041               ikeyecal<EcalIndex.size(); ikeyecal++){
01042             // EcalIndex can have the same cluster save twice (because of the late brem cluster).
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           // End copute TotEnergy
01056 
01057 
01058           // Find extra cluster from e/g importing
01059           for(unsigned int ikeyecal = 0; 
01060               ikeyecal<EcalIndex.size(); ikeyecal++){
01061             // EcalIndex can have the same cluster save twice (because of the late brem cluster).
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                 // check for both KF and GSF
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                     //GSF
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                     //KF
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                     // The new cluster improve both the E/pin?
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         } // END ISOLATION IF
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       // insert the secondary KF tracks
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       // The EcalMap
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             // Check that this cluster is not closer to another ECAL cluster
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             // Check that this cluster is not closer to another ECAL cluster
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           // The first cluster is that one coming from the Gsf. 
01221           // Only for this one is found the HCAL cluster using the Track-HCAL link
01222           // and not the Ecal-Hcal not well tested yet.
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             //  if(localactive[index] == true) {
01232 
01233               // Check that this cluster is not closer to another GSF
01234               // remove in high energetic jets this is dangerous. 
01235 //            std::multimap<double, unsigned> sortedGsf;
01236 //            block.associatedElements( index,  linkData,
01237 //                                      sortedGsf,
01238 //                                      reco::PFBlockElement::GSF,
01239 //                                      reco::PFBlock::LINKTEST_ALL );
01240 //            unsigned jGsf = sortedGsf.begin()->second;
01241 //            if ( jGsf != gsfIs[iEle]) continue; 
01242 
01243               EcalElemsIndex.push_back(index);
01244               localactive[index] = false;
01245               
01246               // }
01247           }
01248           // if the closest ecal cluster has been link with the KF, check KF - HCAL link
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                 // Check that this cluster is not closer to another KF
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           // Find Other Primary Tracks Associated to the same Gsf Clusters
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               // 19 Mar 2010 do not consider here tracks from gamma conv
01285              //  const reco::PFBlockElementTrack * kfTk =  
01286              //  dynamic_cast<const reco::PFBlockElementTrack*>((&elements[index]));         
01287               // DANIELE ?  It is not need because of the local locking 
01288               //   if(kfTk->isLinkedToDisplacedVertex()) continue;
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                   // Check that this track is not closer to another ECAL cluster
01297                   // Not Sure here I need this step
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     }// end type GSF
01316   } // endis there a gsf track
01317   
01318   // ******************* End Link *****************************
01319 
01320   // 
01321   // Below is only for debugging printout 
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   // EndPrintOut
01398   return IsThereAGoodGSFTrack;
01399 }
01400 // This function get the associatedToGsf and associatedToBrems maps and run the final electronID. 
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       // The RefKf is needed to build pure tracking observables
01462       if(assoele_type == reco::PFBlockElement::TRACK) {
01463         const reco::PFBlockElementTrack * KfTk =  
01464           dynamic_cast<const reco::PFBlockElementTrack*>((&elements[(assogsf_index[ielegsf])]));
01465         // 19 Mar 2010 do not consider here track from gamma conv
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           //      Ene_ecalgsf = pfcalib_.energyEm(*clusterRef,ps1Ene,ps2Ene);     
01503           Ene_ecalgsf = pfcalib_.energyEm(*clusterRef,ps1Ene,ps2Ene,ps1,ps2,applyCrackCorrections_);      
01504           //      std::cout << "Test " << Ene_ecalgsf <<  " PS1 / PS2 " << ps1 << " " << ps2 << std::endl;
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       } // end type Ecal
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             // check if it is a compatible cluster also with the gsf track
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       // here build the new BDT observables
01590       
01591       //  ***** Normalization observables ****
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         // **** Pure tracking observables. 
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         // change GsfEl->Pout().pt() as soon the PoutMode is on the GsfTrack DataFormat
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           // Not used, strange behaviour to be checked. And Kf->OuterPt is 
01617           // in track extra. 
01618           //DPtOverPt_kf = 
01619           //  (RefKF->pt() - RefKF->outerPt())/RefKF->pt();
01620           
01621         }
01622         
01623         // **** Tracker-Ecal-Hcal observables 
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         //      std::cout<< " Inserting in extra " << electronExtra_.size() << std::endl;
01649 
01650         // Put cuts and access the BDT output
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         // add output observables 
01681         BDToutput_[cgsf] = mvaValue;
01682         myExtra.setMVA(mvaValue);
01683         electronExtra_.push_back(myExtra);
01684 
01685         // IMPORTANT Additional conditions
01686         if(mvaValue > mvaEleCut_) {
01687           // Check if the ecal cluster is isolated. 
01688           //
01689           // If there is at least one extra track and H/H+E > 0.05 or SumP(ExtraKf)/EGsf or  
01690           // #Tracks > 3 I leave this job to PFAlgo, otherwise I lock those extra tracks. 
01691           // Note:
01692           // The tracks coming from the 4 step of the iterative tracking are not considered, 
01693           // because they can come from secondary converted brems. 
01694           // They are locked to avoid double counting.
01695           // The lock is done in SetActivate function. 
01696           // All this can improved but it is already a good step. 
01697           
01698           
01699           // Find if the cluster is isolated. 
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                 // 19 Mar 2010 do not consider here tracks from gamma conv
01728                 // This should be not needed because the seleted extra tracks from GammaConv
01729                 // will be locked and can not be associated to the ecal elements 
01730                 //if(kfTk->trackType(reco::PFBlockElement::T_FROM_GAMMACONV)) continue;
01731                 
01732 
01733                 reco::TrackRef trackref =  kfTk->trackRef();
01734                 unsigned int Algo = whichTrackAlgo(trackref);
01735                 // iter0, iter1, iter2, iter3 = Algo < 3
01736                 // algo 4,5,6,7
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                   // Check if these extra tracks are HCAL linked
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             // if the Ecluster ~ Ptrack and the extra tracks are HCAL linked 
01773             // the electron is retained and the kf tracks are not locked
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           // This is a pion:
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           // Reject Crazy E/p values... to be understood in the future how to train a 
01801           // BDT in order to avoid to select this bad electron candidates. 
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           // For not-preselected Gsf Tracks ET > 50 GeV, apply dphi preselection
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   } // End Loop on Map1   
01859   return;
01860 }
01861 // This function get the associatedToGsf and associatedToBrems maps and  
01862 // compute the electron 4-mom and set the pf candidate, for
01863 // the gsf track with a BDTcut > mvaEleCut_
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 //   vector<reco::PFCluster> pfClust_vec(0);
01877 //   pfClust_vec.clear();
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     // They should be reset for each gsf track
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     //    vector<unsigned int> assogsf_index =  associatedToGsf_[igsf].second;
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])); // Daniele
01961         const reco::PFBlockElementTrack * KfTk =  
01962           dynamic_cast<const reco::PFBlockElementTrack*>((&elements[(assogsf_index[ielegsf])]));
01963         // 19 Mar 2010 do not consider here track from gamam conv
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         // Important is the PS clusters are not saved before the ecal one, these
01993         // energy are not correctly assigned 
01994         // For the moment I get only the closest PS clusters: this has to be changed
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         //pfClust_vec.push_back((*clust->clusterRef()));
02025 
02026         // The electron RAW energy is the energy of the corrected GSF cluster   
02027         double ps1,ps2;
02028         ps1=ps2=0.;
02029         //      float EE=pfcalib_.energyEm(cl,ps1Ene,ps2Ene);
02030         float EE = pfcalib_.energyEm(cl,ps1Ene,ps2Ene,ps1,ps2,applyCrackCorrections_);    
02031         //      float RawEE = cl.energy();
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           //      std::cout << " PFElectronAlgo / Seed " << EE << std::endl;
02053           RawEene += EE;
02054         }
02055         
02056         // create a photon/electron candidate
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         //            std::cout << " PFElectronAlgo, adding Brem (1) " << EE << std::endl;
02071         // The Raw Ecal energy will be the energy of the basic cluster. 
02072         // It will be the corrected energy without the preshower
02073         cluster_Candidate.setRawEcalEnergy(EE-ps1-ps2);
02074         cluster_Candidate.setPositionAtECALEntrance(math::XYZPointF(cl.position()));
02075         cluster_Candidate.addElementInBlock(blockRef,assogsf_index[ielegsf]);
02076         // store the photon candidate
02077         std::map<unsigned int,std::vector<reco::PFCandidate> >::iterator itcheck=
02078           electronConstituents_.find(cgsf);
02079         if(itcheck==electronConstituents_.end())
02080           {               
02081             // beurk
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       // Important: Add energy from the brems
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             // brem emission is from the considered gsf track
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               //pfClust_vec.push_back(*clusterRef);
02137               // to get a calibrated PS energy 
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               // float RawEE  = clusterRef->energy();
02143               float ceta = clusterRef->position().eta();
02144               // float cphi = clusterRef->position().phi();
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               // Removed 4 March 2009. Florian. The Raw energy is the (corrected) one of the GSF cluster only
02156               //              RawEene += RawEE;
02157               dene+=dE*dE;
02158 
02159               // create a PFCandidate out of it. Watch out, it is for the e/gamma and tau only
02160               // not to be used by the PFAlgo
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               //              std::cout << " PFElectronAlgo, adding Brem " << EE << std::endl;
02174               // yes, EE, we want the raw ecal energy of the daugther to have the same definition
02175               // as the GSF cluster
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               // store the photon candidate
02181               std::map<unsigned int,std::vector<reco::PFCandidate> >::iterator itcheck=
02182                 electronConstituents_.find(cgsf);
02183               if(itcheck==electronConstituents_.end())
02184                 {                 
02185                   // beurk
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     } // End Loop On element associated to the GSF tracks
02200     if (has_gsf) {
02201       
02202       // SuperCluster energy corrections
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       // The calibrations are provided till ET = 200 GeV
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       // charge determination with the majority method
02231       // if the kf track exists: 2 among 3 of supercluster barycenter position
02232       // gsf track and kf track
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       // Think about this... 
02285       if ((nhit_gsf<8) && (has_kf)){
02286         
02287         // Use Hene if some condition.... 
02288         
02289         momentum=momentum_kf;
02290         float Fe=Eene;
02291         float scale= Fe/momentum.E();
02292         
02293         // Daniele Changed
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         // the pf candidate are created: we need to set something more? 
02349         // IMPORTANT -> We need the gsftrackRef, not only the TrackRef??
02350 
02351         if( DebugIDCandidates )
02352           cout << "SetCandidates:: I am before doing candidate " <<endl;
02353         
02354         //vector with the cluster energies (for the extra)
02355         std::vector<float> clusterEnergyVec;
02356         clusterEnergyVec.push_back(RawEene);
02357         clusterEnergyVec.insert(clusterEnergyVec.end(),bremEnergyVec.begin(),bremEnergyVec.end());
02358 
02359         // add the information in the extra
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         // Note the Hcal energy is set but the element is never locked 
02379         temp_Candidate.setHcalEnergy(Hene);  
02380         temp_Candidate.setPs1Energy(ps1TotEne);
02381         temp_Candidate.setPs2Energy(ps2TotEne);
02382         temp_Candidate.setTrackRef(RefKF);   
02383         // This reference could be NULL it is needed a protection? 
02384         temp_Candidate.setGsfTrackRef(RefGSF);
02385         temp_Candidate.setPositionAtECALEntrance(posGsfEcalEntrance);
02386         // Add Vertex
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         // now the special candidate for e/gamma & taus 
02401         reco::PFCandidate extendedElCandidate(temp_Candidate);
02402         // now add the photons to this candidate
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             //      std::cout << " PFElectronAlgo " << nclus << " daugthers to add" << std::endl;
02410             for(unsigned iclus=0;iclus<nclus;++iclus)
02411               {
02412                 extendedElCandidate.addDaughter(theClusters[iclus]);
02413               }
02414           }
02415 //      else
02416 //        {
02417 //          std::cout << " Did not find any photon connected to " <<cgsf << std::endl;
02418 //        }
02419         
02420         allElCandidate_.push_back(extendedElCandidate);
02421         
02422       }
02423       else {
02424         BDToutput_[cgsf] = -1.;   // if the momentum is < 0.5 ID = false, but not sure
02425         // it could be misleading. 
02426         if( DebugIDCandidates ) 
02427           cout << "SetCandidates:: No Candidate Produced because of Pt cut: 0.5 " <<endl;
02428       }
02429     } 
02430     else {
02431       BDToutput_[cgsf] = -1.;  // if gsf ref does not exist
02432       if( DebugIDCandidates ) 
02433         cout << "SetCandidates:: No Candidate Produced because of No GSF Track Ref " <<endl;
02434     }
02435   } // End Loop On GSF tracks
02436   return;
02437 }
02438 // // This function get the associatedToGsf and associatedToBrems maps and for each pfelectron candidate 
02439 // // the element used are set active == false. 
02440 // // note: the HCAL cluster are not used for the moments and also are not locked. 
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     // lock only the elements that pass the BDT cut
02456     if(BDToutput_[cgsf] < mvaEleCut_) continue;
02457 
02458     unsigned int gsf_index =  igsf->first;
02459     active[gsf_index] = false;  // lock the gsf
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       // lock the elements associated to the gsf: ECAL, Brems
02464       active[(assogsf_index[ielegsf])] = false;  
02465       if (assoele_type == reco::PFBlockElement::ECAL) {
02466         unsigned int keyecalgsf = assogsf_index[ielegsf];
02467 
02468         // added protection against fifth step
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         // added locking for conv gsf tracks and kf tracks
02478         if(convGsfTrack_.size() > 0) {
02479           for(unsigned int iconv = 0; iconv < convGsfTrack_.size(); iconv++) {
02480             if(convGsfTrack_[iconv].first == keyecalgsf) {
02481               // lock the GSF track
02482               active[(convGsfTrack_[iconv].second)] = false;
02483               // lock also the KF track associated
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           // lock the elements associated to ECAL: PS1,PS2, for the moment not HCAL
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       } // End if ECAL
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             // lock the ecal cluster associated to the brem
02519             active[(assobrem_index[ibrem])] = false;
02520 
02521             // add protection against fifth step
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             // lock the elements associated to ECAL: PS1,PS2, for the moment not HCAL
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       } // End if BREM    
02541     } // End loop on elements from gsf track
02542   }  // End loop on gsf track  
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 }