---

PFElectronAlgo.cc

Go to the documentation of this file.

//
// Original Author: Daniele Benedetti: daniele.benedetti@cern.ch
//

#include "RecoParticleFlow/PFAlgo/interface/PFElectronAlgo.h"
#include "DataFormats/ParticleFlowReco/interface/PFBlockElementGsfTrack.h"
#include "DataFormats/ParticleFlowReco/interface/PFBlockElementCluster.h"
#include "DataFormats/ParticleFlowReco/interface/PFBlockElementBrem.h"
#include "DataFormats/ParticleFlowReco/interface/PFRecHitFraction.h"
#include "DataFormats/ParticleFlowReco/interface/PFRecHitFwd.h" 
#include "DataFormats/ParticleFlowReco/interface/PFRecHit.h"
#include "DataFormats/ParticleFlowReco/interface/PFCluster.h"
#include "DataFormats/GsfTrackReco/interface/GsfTrack.h"
#include "DataFormats/GsfTrackReco/interface/GsfTrackFwd.h"
#include "RecoParticleFlow/PFClusterTools/interface/PFEnergyCalibration.h"
#include "RecoParticleFlow/PFClusterTools/interface/PFEnergyResolution.h"

using namespace std;
using namespace reco;
PFElectronAlgo::PFElectronAlgo(const std::vector<double> chi2values,
                               const double mvaEleCut,
                               string mvaWeightFileEleID):
  chi2values_(chi2values),
  mvaEleCut_(mvaEleCut)
{
  // Set the tmva reader
  tmvaReader_ = new TMVA::Reader();
  tmvaReader_->AddVariable("TEtotPinMode",&TEtotPinMode);
  tmvaReader_->AddVariable("TEGsfPoutMode",&TEGsfPoutMode);
  tmvaReader_->AddVariable("TDiffEtaGSFEcalEle",&TDiffEtaGSFEcalEle);
  tmvaReader_->AddVariable("log(TSigmaEtaEtaEle)",&TSigmaEtaEtaEle);
  tmvaReader_->AddVariable("TChiEcalEle",&TChiEcalEle);
  tmvaReader_->AddVariable("TNumBrem",&TNumBrem);
  tmvaReader_->AddVariable("TSumEBremPinPoutMode",&TSumEBremPinPoutMode);
  tmvaReader_->AddVariable("TLateBrem",&TLateBrem);
  tmvaReader_->AddVariable("TEarlyBrem",&TEarlyBrem);
  tmvaReader_->BookMVA("BDT",mvaWeightFileEleID.c_str());
}
void PFElectronAlgo::RunPFElectron(const reco::PFBlockRef&  blockRef,
                                   std::vector<bool>& active){
  // General note: all the cout will be removed. 

  // the maps are initialized 
  AssMap associatedToGSF(0);
  AssMap associatedToBrems(0);
  // this function fill the previous maps: gsf and associated elements
  // and brem and associated elements
  // To be noted: for the moment the function FindClosestElement is based
  // only on geometrical link (by chi2) but Jonathan is working on a 
  // geometrical + energetic link, see presentation:
  // https://indico.cern.ch/getFile.py/access?contribId=0&resId=1&materialId=slides&confId=33548

  
  GsfTrackSingleEcal_.clear();
  GsfTrackSingleEcal_ = active;
  MultipleGsfTrackAssociation(blockRef);

  bool blockHasGSF = SetLinks(blockRef,associatedToGSF,associatedToBrems,active);
  
  if (blockHasGSF) {
    // If at least a gsf pre-id track has been connected to an ecal cluster
    // is run the electron - ID
    // When will have also pure tracking-ID this will be removed. 
    
    BDToutput_.clear();
    // For each GSF track is initialized a BDT value = -1
    BDToutput_.assign(associatedToGSF.size(),-1.);
    
    // The FinalID is run and BDToutput is internally changed. 
    SetIDOutputs(blockRef,associatedToGSF,associatedToBrems);  
    elCandidate_.clear();
    // For each GSF track that pass the BDT configurable cut is run the setcandidates.
    // This function computed the initial electron 4-momentum and fill the vector 
    // of pf candidates
    SetCandidates(blockRef,associatedToGSF,associatedToBrems);
    if (elCandidate_.size() > 0 ){
      isvalid_ = true;
      // if there is at least a candidate the active vector is modified
      // setting = false all the elements used to build the candidate
      SetActive(blockRef,associatedToGSF,associatedToBrems,active);
      }
  } // endif blockHasGSF
}
// For the moment this function just get the closest element by chi2: 
// in future the link could be geometrical + energetic so it is good to have
// this separated function
uint PFElectronAlgo::FindClosestElement(const uint iele,
                                        std::multimap<double, unsigned>& Elems, 
                                        float& chi2cut,
                                        std::vector<bool>& active,
                                        const reco::PFBlockRef&  blockRef) {
  
  const reco::PFBlock& block = *blockRef;
  PFBlock::LinkData linkData =  block.linkData();     
  uint FoundIndex = 100000;
  if (Elems.size() > 0) {
    std::multimap<double, unsigned>::iterator ie = Elems.begin();
    // float chi2  = ie->first;
    uint index = ie->second;
    float chi2 =  block.chi2(iele,index,linkData);
    if (chi2 < chi2cut && active[index] == true) {
      chi2cut = chi2;
      FoundIndex = index;
    }
  }
  return FoundIndex;
}
// Like FindClosestElement, but check if the cluster has been alredy associated 
// to the gsf track. It needs to understand if there is a late brems
uint PFElectronAlgo::FindClosestElementCond(const uint iele,
                                            std::multimap<double, unsigned>& Elems, 
                                            float& chi2cut,float chi2Comp, 
                                            uint indexComp,
                                            std::vector<bool>& active,
                                            const reco::PFBlockRef&  blockRef) {
  const reco::PFBlock& block = *blockRef;
  PFBlock::LinkData linkData =  block.linkData();     
  uint FoundIndex = 100000;
  for(std::multimap<double, unsigned>::iterator ie = Elems.begin(); ie != Elems.end(); ++ie ) {
    if (ie->second == indexComp && (ie->first + 5.) > chi2Comp) continue; 
    //    double   chi2  = ie->first;
    unsigned index = ie->second;
    float chi2 =  block.chi2(iele,index,linkData);
    if (chi2 < chi2cut && active[index] == true ) {
      chi2cut = chi2;
      FoundIndex = index;
    }
    break;
  }
  return FoundIndex;
}
// Solve the double counting: for each brem the closest ecal cluster is found.
// If an ecal cluster has been associated to several brems the closest is chosen. 
bool PFElectronAlgo::findLinkPairs(const vector< pair<unsigned, double> >& MinRow , 
                                   const vector<unsigned>& brem_it, 
                                   const vector<unsigned>& clust_it, 
                                   map<unsigned, unsigned>& finallink, 
                                   float chi2cut){


  if (MinRow.size() >0 ) {

    for (unsigned i1=0; i1< clust_it.size(); i1++) {
      double final_chi2 = chi2cut;  
      unsigned final_brem_index = 100000;
      unsigned final_ecal_index = 100000;      
      for(unsigned i2=0;i2<MinRow.size();i2++) {
        if (MinRow[i2].first == clust_it[i1]) {
          if (MinRow[i2].second < final_chi2) {
            final_chi2 = MinRow[i2].second;
            final_brem_index = brem_it[i2];
            final_ecal_index = clust_it[i1];
          }
        }
      }
      if (final_chi2 < chi2cut) {  
        finallink.insert(make_pair(final_brem_index,final_ecal_index));
      }
    }
    if (finallink.size() > 0)  return true;
    else return false;
  }
  else {
    return false;
  }
}
// The maps associatedToGSF and associatedToBrems are filled. 
void PFElectronAlgo::MultipleGsfTrackAssociation(const reco::PFBlockRef&  blockRef){
  

  const reco::PFBlock& block = *blockRef;
  const edm::OwnVector< reco::PFBlockElement >&  elements = block.elements();
  PFBlock::LinkData linkData =  block.linkData(); 
  for(unsigned iEle=0; iEle<elements.size(); iEle++) {    
    PFBlockElement::Type type = elements[iEle].type();
    if (type == reco::PFBlockElement::ECAL ) {
      std::multimap<double, unsigned> gsfTrack;
      block.associatedElements( iEle,  linkData,
                                gsfTrack ,
                                reco::PFBlockElement::GSF,
                                reco::PFBlock::LINKTEST_ALL );
      if (gsfTrack.size() > 1.) {
        for(std::multimap<double, unsigned>::iterator ie = gsfTrack.begin(); 
            ie != gsfTrack.end(); ++ie ) {
          unsigned index_gsf = ie->second;
          GsfTrackSingleEcal_[index_gsf] = false;
        }
      }
    }
  }
  return;
}
// The maps associatedToGSF and associatedToBrems are filled. 
bool PFElectronAlgo::SetLinks(const reco::PFBlockRef&  blockRef,
                              AssMap& associatedToGSF_,
                              AssMap& associatedToBrems_,
                              std::vector<bool>& active){
  uint CutIndex = 100000;
  double CutGSFECAL = chi2values_[0] ;   // Standard Value should be 900 
  double CutBremECAL = chi2values_[1];   // Standard Value should be 25
  double CutGSFHCAL = chi2values_[2];    // Standard Value sould be 100, to be studied
  double CutBremHCAL = chi2values_[3];   // Standard Value sould be 25, to be studied
  double CutGSFPS = chi2values_[4];      // Standard Value sould be 100, to be studied
  double CutBremPS = chi2values_[5];     // Standard Value sould be 25, to be studied
  bool DebugSetLinks = false;
  
  const reco::PFBlock& block = *blockRef;
  const edm::OwnVector< reco::PFBlockElement >&  elements = block.elements();
  PFBlock::LinkData linkData =  block.linkData();  
  
  bool IsThereAGSFTrack = false;
  for(unsigned iEle=0; iEle<elements.size(); iEle++) {
     
     std::map<unsigned, unsigned> FinalBremECALIndex;
     std::map<unsigned, unsigned> FinalBremHCALIndex; 
     std::map<unsigned, unsigned> FinalBremPS1Index;  
     std::map<unsigned, unsigned> FinalBremPS2Index;  
     unsigned ECALGSF_index = CutIndex;
     unsigned GSF_index = CutIndex;
     unsigned KFGSF_index = CutIndex;
     unsigned HCALGSF_index = CutIndex;
     unsigned PS1GSF_index = CutIndex;
     unsigned PS2GSF_index = CutIndex;

     PFBlockElement::Type type = elements[iEle].type();
     if (type == reco::PFBlockElement::GSF ) {
       IsThereAGSFTrack = true;
       GSF_index = iEle;

       if (!active[iEle]) continue;  // not use if locked by muons
       if (DebugSetLinks) cout << "PFElectron:: The Block " << block << endl;
       // Stop if more GSF tracks are associated to the same clusters.
       // to be modified in the future. 
       if(!GsfTrackSingleEcal_[iEle]) continue;
       

       // Search for elements associated to the GSF track
       // Find Associated KF Tracks
       std::multimap<double, unsigned> kfTrack;
       block.associatedElements( iEle,  linkData,
                                 kfTrack ,
                                 reco::PFBlockElement::TRACK,
                                 reco::PFBlock::LINKTEST_ALL );
       if (kfTrack.size() > 0) {
         multimap<double, unsigned>::iterator itkf = kfTrack.begin();
         KFGSF_index = itkf->second;
       }


       // Find Associated ECAL Clusters
       std::multimap<double, unsigned> ecalElems;
       block.associatedElements( iEle,  linkData,
                                 ecalElems ,
                                 reco::PFBlockElement::ECAL,
                                 reco::PFBlock::LINKTEST_ALL );     
       float chi2ECALGSFFinal = CutGSFECAL;   
       ECALGSF_index = FindClosestElement(iEle,ecalElems,
                                          chi2ECALGSFFinal,
                                          active,blockRef);
   
       if (DebugSetLinks) {
         if (ecalElems.size() > 0) {
           cout << "PFElectron:: ECALGSF_index " << ECALGSF_index  << " chi2 " << chi2ECALGSFFinal << endl;
         }
         else {
           cout << "PFElectron:: No Link! " << endl;
         }
       }

       // Find Associated HCAL Clusters
       std::multimap<double, unsigned> hcalElems;
       block.associatedElements( iEle,  linkData,
                                 hcalElems ,
                                 reco::PFBlockElement::HCAL,
                                 reco::PFBlock::LINKTEST_ALL );     
       float chi2HCALGSFFinal = CutGSFHCAL;
       HCALGSF_index = FindClosestElement(iEle,hcalElems,
                                          chi2HCALGSFFinal,
                                          active,blockRef);
       
       
       // Find Associated PS1 Clusters
       std::multimap<double, unsigned> ps1Elems;
       block.associatedElements( iEle,  linkData,
                                 ps1Elems ,
                                 reco::PFBlockElement::PS1,
                                 reco::PFBlock::LINKTEST_ALL );     
       float chi2PS1GSFFinal = CutGSFPS;
       PS1GSF_index = FindClosestElement(iEle,ps1Elems,
                                         chi2PS1GSFFinal,
                                         active,blockRef);
       
       
       // Find Associated PS2 Clusters
       std::multimap<double, unsigned> ps2Elems;
       block.associatedElements( iEle,  linkData,
                                 ps2Elems ,
                                 reco::PFBlockElement::PS1,
                                 reco::PFBlock::LINKTEST_ALL );     
       float chi2PS2GSFFinal = CutGSFPS;
       PS2GSF_index = FindClosestElement(iEle,ps2Elems,
                                         chi2PS2GSFFinal,
                                         active,blockRef); 
       // End First Part Search for elements associated to the GSF track
       
       
       // Second Part search for brems and elements associated to the brems 
       std::multimap<double, unsigned> Brems;
       
       block.associatedElements( iEle,  linkData,
                                 Brems ,
                                 reco::PFBlockElement::BREM,
                                 reco::PFBlock::LINKTEST_ALL );
       
       

       vector<unsigned> brems_index(0);
       vector<unsigned> mainecal_index(0);
       vector< pair<unsigned, double> >  BremEcal(0);
 
       for(std::multimap<double, unsigned>::iterator ieb = Brems.begin(); 
           ieb != Brems.end(); ++ieb ) {
         
         
         // ECAL associated to the brems
         unsigned index_brem = ieb->second;
         if (DebugSetLinks) cout << "PFElectron:: index brem " << index_brem  << endl; 

         std::multimap<double, unsigned> ecalBremsElems;
         block.associatedElements( index_brem,  linkData,
                                   ecalBremsElems,
                                   reco::PFBlockElement::ECAL,
                                   reco::PFBlock::LINKTEST_ALL );
         float chi2EcalBremFinal = CutBremECAL;
         unsigned ECALBREM_index = FindClosestElementCond(index_brem,ecalBremsElems,chi2EcalBremFinal,
                                                          chi2ECALGSFFinal,ECALGSF_index,
                                                          active,blockRef);
         
         if (ECALBREM_index < CutIndex){          
           bool NonAlreadyAss = true;
           if (BremEcal.size() > 0) {
             for (uint ieb = 0; ieb < BremEcal.size(); ieb++){
               if (ECALBREM_index == BremEcal[ieb].first) {
                 NonAlreadyAss = false;
                 break;
               }
             }
           }
           BremEcal.push_back(make_pair(ECALBREM_index,chi2EcalBremFinal));
           brems_index.push_back(index_brem);
           if (NonAlreadyAss) mainecal_index.push_back(ECALBREM_index);
         }
       }
       
       // This solves the possible double counting for ecal clusters associated to brems. 
       bool foundEcalLink = findLinkPairs(BremEcal,brems_index,mainecal_index, 
                                          FinalBremECALIndex,CutBremECAL);
       
       
       vector<unsigned> element_brems(0);
       if (foundEcalLink) {
         for (map<unsigned, unsigned>::iterator it = FinalBremECALIndex.begin(); 
              it != FinalBremECALIndex.end(); ++it ) {
           unsigned index_bremassoecal = it->first;
           element_brems.push_back(index_bremassoecal);
           
            if (DebugSetLinks) cout <<  "PFElectron:: linked index brem " << it->first  
                                    <<" index ecal " << it->second << endl; 

           // HCAL associated to the brems
           std::multimap<double, unsigned> hcalBrems;
           block.associatedElements( index_bremassoecal,  linkData,
                                     hcalBrems ,
                                     reco::PFBlockElement::HCAL,
                                     reco::PFBlock::LINKTEST_ALL );
           float chi2HcalBremFinal = CutBremHCAL;
           unsigned HCALBREM_index = FindClosestElementCond(index_bremassoecal,hcalBrems,chi2HcalBremFinal,
                                                            chi2HCALGSFFinal,HCALGSF_index,
                                                            active,blockRef);
           if (chi2HcalBremFinal < CutBremHCAL) 
             FinalBremHCALIndex.insert(make_pair(index_bremassoecal,HCALBREM_index));
           
           

           // PS1 associated to the brems
           std::multimap<double, unsigned> PS1Brems;
           block.associatedElements( index_bremassoecal,  linkData,
                                     PS1Brems ,
                                     reco::PFBlockElement::PS1,
                                     reco::PFBlock::LINKTEST_ALL );
           float chi2PS1BremFinal = CutBremPS;
           unsigned PS1BREM_index  = FindClosestElementCond(index_bremassoecal,PS1Brems,chi2PS1BremFinal,
                                                            chi2PS1GSFFinal,PS1GSF_index,
                                                            active,blockRef);
           if (chi2PS1BremFinal < CutBremPS) 
             FinalBremPS1Index.insert(make_pair(index_bremassoecal,PS1BREM_index));
           
           

           // PS2 associated to the brems
           std::multimap<double, unsigned> PS2Brems;
           block.associatedElements( index_bremassoecal,  linkData,
                                     PS2Brems ,
                                     reco::PFBlockElement::PS2,
                                     reco::PFBlock::LINKTEST_ALL );
           float chi2PS2BremFinal = CutBremPS;
           unsigned PS2BREM_index = FindClosestElementCond(index_bremassoecal,PS2Brems,chi2PS2BremFinal,
                                                           chi2PS2GSFFinal,PS2GSF_index,
                                                           active,blockRef);
           if (chi2PS2BremFinal < CutBremPS) 
             FinalBremPS2Index.insert(make_pair(index_bremassoecal,PS2BREM_index));
         }
       }
       
            
       // Save an AssMap with elements associated to the brem track
       // Note the order is important in setCandidate and setIDOutput 
       vector<unsigned> assoel(0);
       if (KFGSF_index < CutIndex) assoel.push_back(KFGSF_index);
       if (PS1GSF_index  < CutIndex) assoel.push_back(PS1GSF_index);
       if (PS2GSF_index  < CutIndex) assoel.push_back(PS2GSF_index);
       if (ECALGSF_index < CutIndex)  assoel.push_back(ECALGSF_index); 
       if (HCALGSF_index < CutIndex) assoel.push_back(HCALGSF_index);
  
       for (map<unsigned, unsigned>::iterator it = FinalBremECALIndex.begin(); 
            it != FinalBremECALIndex.end(); ++it ) {
         unsigned index_brem = it->first;
         assoel.push_back(index_brem);
       }
       associatedToGSF_.push_back(make_pair(GSF_index,assoel));
       
       // Save an AssMap with elements associated to the brem track
       // Note the order is important in setCandidate and setIDOutput 
       bool DebugBrem = false;
       for(unsigned iBrem=0; iBrem<element_brems.size(); iBrem++) {
         vector<unsigned> assobrem(0);
         for (map<unsigned, unsigned>::iterator it = FinalBremPS1Index.begin(); 
              it != FinalBremPS1Index.end(); ++it ) {
           if (it->first == element_brems[iBrem]) {
             unsigned index_clus = it->second;
             assobrem.push_back(index_clus);
             if(DebugBrem) cout << " In associate BREM " <<  element_brems[iBrem] 
                                << " PS1 Cluster Associated to BREM " 
                                <<  index_clus << endl;
           }
         }
         for (map<unsigned, unsigned>::iterator it = FinalBremPS2Index.begin(); 
              it != FinalBremPS2Index.end(); ++it ) {
           if (it->first == element_brems[iBrem]) {
             unsigned index_clus = it->second;
             assobrem.push_back(index_clus);
             if(DebugBrem) cout << " In associate BREM " <<  element_brems[iBrem] 
                                << " PS2 Cluster Associated to BREM " 
                                <<  index_clus << endl;
           }
         }

         for (map<unsigned, unsigned>::iterator it = FinalBremECALIndex.begin(); 
              it != FinalBremECALIndex.end(); ++it ) {
           if (it->first == element_brems[iBrem]) {
             unsigned index_clus = it->second;
             assobrem.push_back(index_clus);
             // if(ev == 1957) 
             if(DebugBrem) cout << " In associate BREM " <<  element_brems[iBrem] 
                                << " ECAL Cluster Associated to BREM " 
                                <<  index_clus << endl;
           }
         }
         for (map<unsigned, unsigned>::iterator it = FinalBremHCALIndex.begin(); 
              it != FinalBremHCALIndex.end(); ++it ) {
           if (it->first == element_brems[iBrem]) {
             unsigned index_clus = it->second;
             assobrem.push_back(index_clus);
             if(DebugBrem) cout << " In associate BREM " <<  element_brems[iBrem] 
                                << " HCAL Cluster Associated to BREM " 
                                <<  index_clus << endl;
           }
         }
         associatedToBrems_.push_back(make_pair(element_brems[iBrem],assobrem));
       } // End Loop On Brems for GSF track               
     }// end type GSF
  }// End Loop on Block's element

  return IsThereAGSFTrack;
}
// This function get the associatedToGSF and associatedToBrems maps and run the final electronID. 
void PFElectronAlgo::SetIDOutputs(const reco::PFBlockRef&  blockRef,
                                  const AssMap& associatedToGSF_,
                                  const AssMap& associatedToBrems_){
  PFEnergyCalibration pfcalib_;
  const reco::PFBlock& block = *blockRef;
  PFBlock::LinkData linkData =  block.linkData();     
  const edm::OwnVector< reco::PFBlockElement >&  elements = block.elements();
 
  for (uint igsf = 0; igsf<associatedToGSF_.size();igsf++) {
    uint index_gsf =  associatedToGSF_[igsf].first;
    PFBlockElement::Type typegsf = elements[index_gsf].type();

    // Find position in the ECAl surface

    double GSF_Ein   = 0.;
    double GSF_Eout  = 0.;
    double GSF_etaout = 0.;

    bool DebugIDOutputs = false;

    if (typegsf == reco::PFBlockElement::GSF) {  // Just to be sure, but not needed
      const reco::PFBlockElementGsfTrack * GsfEl  =  
        dynamic_cast<const reco::PFBlockElementGsfTrack*>((&elements[index_gsf]));
      math::XYZTLorentzVector p_out = GsfEl->Pout();  // can be mode or mean: see config file xyz

      GSF_Eout = p_out.t();
    
      reco::GsfTrackRef RefGSF = GsfEl->GsftrackRef();
      if (RefGSF.isNonnull()) {
        float m_el=0.00051;
        GSF_Ein = sqrt(RefGSF->pMode()*
                       RefGSF->pMode()+m_el*m_el);
      }

      // Gsfouter position on the ECAL surface
      const math::XYZPointF& PosEcal = GsfEl->positionAtECALEntrance();
      GSF_etaout = PosEcal.eta();
 
      if (DebugIDOutputs)  cout << " setIdOutput! GSF Track " << GSF_Ein  
                                << " ECAL eta,phi " << GSF_etaout   
                                <<", " << PosEcal.phi() << endl;
    } // End if typegsf == GSF
    
    
    vector<uint> index_assogsf =  associatedToGSF_[igsf].second;

    double GSFClus_energy = 0.;
    double Sigmaetaeta = 0.;
    double DiffEtaGSFClust = 0.;
    unsigned GSFClus_index = 100000;
    double chi2GSFECAL = 0.;
    double Tot_BremClustEne = 0.;
    bool LateBrem = false;
    bool EarlyBrem = false;
    int NumBrem = 0;
    
    // The final Observables are build starting from the gsf and brem maps 

    for  (unsigned ielegsf=0;ielegsf<index_assogsf.size();ielegsf++) {
      PFBlockElement::Type typeassoele = elements[(index_assogsf[ielegsf])].type();
      vector<double> ps1Ene(0);
      vector<double> ps2Ene(0);
      // Important is the PS clusters are not saved before the ecal one, these
      // energy are not correctly assigned 
      // For the moment I get only the closest PS clusters: this has to be changed
      if  (typeassoele == reco::PFBlockElement::PS1) {  
        PFClusterRef  psref = elements[(index_assogsf[ielegsf])].clusterRef();
        ps1Ene.push_back(psref->energy());
      }
      if  (typeassoele == reco::PFBlockElement::PS2) {  
        PFClusterRef  psref = elements[(index_assogsf[ielegsf])].clusterRef();
        ps2Ene.push_back(psref->energy());
      }
      if  (typeassoele == reco::PFBlockElement::ECAL) {
        GSFClus_index = index_assogsf[ielegsf];
        reco::PFClusterRef clusterRef = elements[(index_assogsf[ielegsf])].clusterRef();
        //GSFClus_energy = clusterRef->energy();
        GSFClus_energy = pfcalib_.energyEm(*clusterRef,ps1Ene,ps2Ene);

        chi2GSFECAL = block.chi2(index_gsf,index_assogsf[ielegsf],linkData);
        double GSFClus_eta = clusterRef->position().eta();
        if (DebugIDOutputs)  cout << " setIdOutput! GSF ECAL Cluster E " << GSFClus_energy
                                  << " eta,phi " << GSFClus_eta
                                  <<", " <<  clusterRef->position().phi() << endl;

        DiffEtaGSFClust = GSFClus_eta - GSF_etaout;

        const vector< reco::PFRecHitFraction >& PFRecHits =  clusterRef->recHitFractions();
        double SeedClusEnergy = -1.;
        unsigned SeedDetID = 0;
        double SeedEta = -1.;
        double SeedPhi = -1.;
        for ( vector< reco::PFRecHitFraction >::const_iterator it = PFRecHits.begin(); 
              it != PFRecHits.end(); ++it) {
          const PFRecHitRef& RefPFRecHit = it->recHitRef(); 
          double energyRecHit = RefPFRecHit->energy();
          if (energyRecHit > SeedClusEnergy) {
            SeedClusEnergy = energyRecHit;
            SeedEta = RefPFRecHit->position().eta();
            SeedPhi =  RefPFRecHit->position().phi();
            SeedDetID = RefPFRecHit->detId();
          }
        }
        
        for ( vector< reco::PFRecHitFraction >::const_iterator it = PFRecHits.begin(); 
              it != PFRecHits.end(); ++it) {
          
          const PFRecHitRef& RefPFRecHit = it->recHitRef(); 
          double energyRecHit = RefPFRecHit->energy();
          if (RefPFRecHit->detId() != SeedDetID) {
            float diffEta =  RefPFRecHit->position().eta() - SeedEta;
            Sigmaetaeta += (diffEta*diffEta) * (energyRecHit/SeedClusEnergy);
          }
        }
        if (Sigmaetaeta == 0.) Sigmaetaeta = 0.00000001;
      } // End if on ECAL element
      // Add Preshower info for the GSF-PS-ECAL FinalID

      if  (typeassoele == reco::PFBlockElement::BREM) {
        for (uint ibrem = 0; ibrem < associatedToBrems_.size(); ibrem++){
          vector<unsigned> index_assobrem = (associatedToBrems_[ibrem]).second;
          unsigned index_brem = (associatedToBrems_[ibrem]).first;        
          if (index_brem == index_assogsf[ielegsf]) {   // Note: important: check that 
            // brem emission is from the considered gsf track
            const reco::PFBlockElementBrem * BremEl  =  
              dynamic_cast<const reco::PFBlockElementBrem*>((&elements[index_brem]));
            unsigned BremTrajP = BremEl->indTrajPoint();
            unsigned TrajPos = BremTrajP-2;
//          double BremDP = (BremEl->DeltaP())*(-1.); // not used yet
//          double BremSDP = BremEl->SigmaDeltaP(); // not used yet
            if (TrajPos <= 3) EarlyBrem = true;
            for (unsigned ielebrem=0;ielebrem<index_assobrem.size();ielebrem++) {
              vector<double> ps1EneFromBrem(0);
              vector<double> ps2EneFromBrem(0);
              // Important is the PS clusters are not saved before the ecal one, these
              // energy are not correctly assigned 
              // For the moment I get only the closest PS clusters: this has to be changed
              if (elements[(index_assobrem[ielebrem])].type() == reco::PFBlockElement::PS1) {
                PFClusterRef  psref = elements[(index_assobrem[ielebrem])].clusterRef();
                ps1EneFromBrem.push_back(psref->energy());
              }
              if (elements[(index_assobrem[ielebrem])].type() == reco::PFBlockElement::PS2) {
                PFClusterRef  psref = elements[(index_assobrem[ielebrem])].clusterRef();
                ps2EneFromBrem.push_back(psref->energy());
              }        
              if (elements[(index_assobrem[ielebrem])].type() == reco::PFBlockElement::ECAL) {
                // check if it is a compatible cluster also with the gsf track
                if( index_assobrem[ielebrem] !=  GSFClus_index) {
                  reco::PFClusterRef clusterRef = 
                    elements[(index_assobrem[ielebrem])].clusterRef();
                  //double BremClus_energy = clusterRef->energy();
                  double BremClus_energy = 
                    pfcalib_.energyEm(*clusterRef,ps1EneFromBrem,ps2EneFromBrem);
                  double chi2BREMECAL = 
                    block.chi2(index_brem,(index_assobrem[ielebrem]),linkData);
                   if (DebugIDOutputs) cout << " setIdOutput:: GSF BREM Cluster " 
                                            <<  BremClus_energy << " eta,phi " 
                                            <<  clusterRef->position().eta()  <<", " 
                                            << clusterRef->position().phi() 
                                            << " chi " << sqrt(chi2BREMECAL) << endl;
                  Tot_BremClustEne += BremClus_energy;
                  NumBrem++;
                }
                else LateBrem = true;
              }
            }
          }       
        }
      } // End if on BREM elements
    } // End Loop on GSF elements
    if (GSFClus_energy > 0.) {
      TEtotPinMode        =  (GSFClus_energy + Tot_BremClustEne) / GSF_Ein;  
      if (TEtotPinMode > 5.) TEtotPinMode = 5.;
      TEGsfPoutMode       =  GSFClus_energy/GSF_Eout;
      if (TEGsfPoutMode > 5.) TEGsfPoutMode = 5.;
      TDiffEtaGSFEcalEle  = fabs(DiffEtaGSFClust);
      if (TDiffEtaGSFEcalEle > 0.2)  TDiffEtaGSFEcalEle = 0.2;
      TSigmaEtaEtaEle     = Sigmaetaeta;
      if (TSigmaEtaEtaEle > 0.05) TSigmaEtaEtaEle = 0.05;
      TSigmaEtaEtaEle = log(TSigmaEtaEtaEle);
      
      TChiEcalEle         = sqrt(chi2GSFECAL); 
      
      TNumBrem = NumBrem;
      TSumEBremPinPoutMode = 0.;
      TLateBrem = -1;
      TEarlyBrem = -1;
      if (TNumBrem > 0) {
        TSumEBremPinPoutMode = Tot_BremClustEne /(GSF_Ein - GSF_Eout);
        if (TSumEBremPinPoutMode < 0.) TSumEBremPinPoutMode = 0.;
        if (TSumEBremPinPoutMode > 5.) TSumEBremPinPoutMode = 5.;
        
        if (LateBrem == true) TLateBrem = 1;
        else TLateBrem = 0;
        if (EarlyBrem == true) TEarlyBrem = 1;
        else TEarlyBrem = 0;
      }

      if (DebugIDOutputs) cout << " variables " << 
        " TEtotPinMode "  << TEtotPinMode << "  " << 
        " TEGsfPoutMode "  << TEGsfPoutMode << "  " << 
        " TDiffEtaGSFEcalEle "  << TDiffEtaGSFEcalEle  << "  " << 
        " TSigmaEtaEtaEle "  << TSigmaEtaEtaEle << "  " << 
        " TChiEcalEle "  << TChiEcalEle << "  " << 
        " TNumBrem "  << TNumBrem << "  " << 
        " TSumEBremPinPoutMode "  << TSumEBremPinPoutMode << "  " << 
        " TLateBrem "  << TLateBrem << "  " << 
        " TEarlyBrem "  << TEarlyBrem  << "  " << endl;
      
      // the BDT is run 
      double mvaValue = tmvaReader_->EvaluateMVA("BDT");
      //Double_t mvaValue  = 0.5;

      // This condition has to be modified as soon we have
      // better BDT observables and training. 
      if (TEtotPinMode > 1.5 && TEarlyBrem != 1) {BDToutput_[igsf] = -1.;}
      else {BDToutput_[igsf] = mvaValue;}
  
      
      if (DebugIDOutputs) cout << "setIdOutput:: BDT output " 
                               << igsf << " mvaValue " 
                               << BDToutput_[igsf] << endl;
    }
  } // End Loop on Map1 

  return;
}
// This function get the associatedToGSF and associatedToBrems maps and  
// compute the electron 4-mom and set the pf candidate, for
// the gsf track with a BDTcut > mvaEleCut_
void PFElectronAlgo::SetCandidates(const reco::PFBlockRef&  blockRef,
                                   const AssMap& associatedToGSF_,
                                   const AssMap& associatedToBrems_) {

  const reco::PFBlock& block = *blockRef;
  PFBlock::LinkData linkData =  block.linkData();     
  const edm::OwnVector< reco::PFBlockElement >&  elements = block.elements();
  PFEnergyResolution pfresol_;
  PFEnergyCalibration pfcalib_;
  bool DebugIDCandidates = false;

  for (uint igsf = 0; igsf<associatedToGSF_.size();igsf++) {
    uint index_gsf =  associatedToGSF_[igsf].first;
    if(BDToutput_[igsf] < mvaEleCut_) continue;  
    // Important enter only if good cut: to set configurable

    // They should be reset for each gsf track
    int eecal=0;
    int hcal=0;
    int charge =0; 
    bool goodphi=true;
    math::XYZTLorentzVector momentum_kf,momentum_gsf,momentum,momentum_mean;
    float dpt=0; float dpt_kf=0; float dpt_gsf=0; float dpt_mean=0;
    float Eene=0; float dene=0; float Hene=0.;
    float de_gs, de_me, de_kf; 
    float m_el=0.00051;
    int nhit_kf=0; int nhit_gsf=0;
    bool has_gsf=false;
    bool has_kf=false;
    math::XYZTLorentzVector newmomentum;
    float ps1TotEne = 0;
    float ps2TotEne = 0;
    vector<uint> elementsToAdd(0);
    reco::GsfTrackRef RefGSF;  
    reco::TrackRef RefKF;  


    PFBlockElement::Type typegsf = elements[index_gsf].type();
    if (typegsf == reco::PFBlockElement::GSF) { 
      elementsToAdd.push_back(index_gsf);
      const reco::PFBlockElementGsfTrack * GsfEl  =  
        dynamic_cast<const reco::PFBlockElementGsfTrack*>((&elements[index_gsf]));
      RefGSF = GsfEl->GsftrackRef();
      if (RefGSF.isAvailable()) {
        if (RefGSF.isNonnull()) {

          has_gsf=true;

          charge= RefGSF->chargeMode();
          nhit_gsf= RefGSF->hitPattern().trackerLayersWithMeasurement();

          momentum_gsf.SetPx(RefGSF->pxMode());
          momentum_gsf.SetPy(RefGSF->pyMode());
          momentum_gsf.SetPz(RefGSF->pzMode());
          float ENE=sqrt(RefGSF->pMode()*
                         RefGSF->pMode()+m_el*m_el);

          if( DebugIDCandidates ) 
            cout << "SetCandidates:: GsfTrackRef: Ene " << ENE 
                 << " charge " << charge << endl;

          momentum_gsf.SetE(ENE);       
          dpt_gsf=RefGSF->ptModeError()*
            (RefGSF->pMode()/RefGSF->ptMode());
          
          momentum_mean.SetPx(RefGSF->px());
          momentum_mean.SetPy(RefGSF->py());
          momentum_mean.SetPz(RefGSF->pz());
          float ENEm=sqrt(RefGSF->p()*
                          RefGSF->p()+m_el*m_el);
          momentum_mean.SetE(ENEm);       
          dpt_mean=RefGSF->ptError()*
            (RefGSF->p()/RefGSF->pt());  
        }
        else {
          if( DebugIDCandidates ) 
            cout <<  "SetCandidates:: !!!!  NULL GSF Track Ref " << endl;       
        } 
      }  
      else {
        if( DebugIDCandidates ) 
          cout <<  "SetCandidates:: !!!! NotAvailable GSF Track Ref " << endl;
      }
    }

    vector<uint> index_assogsf =  associatedToGSF_[igsf].second;
    unsigned GSFClus_index = 100000;
    for  (unsigned ielegsf=0;ielegsf<index_assogsf.size();ielegsf++) {
      PFBlockElement::Type typeassoele = elements[(index_assogsf[ielegsf])].type();
      if  (typeassoele == reco::PFBlockElement::TRACK) {
        elementsToAdd.push_back((index_assogsf[ielegsf])); // Daniele
        const reco::PFBlockElementTrack * KfTk =  
          dynamic_cast<const reco::PFBlockElementTrack*>((&elements[(index_assogsf[ielegsf])]));        
        RefKF = KfTk->trackRef();
        if (RefKF.isAvailable()) {
          if (RefKF.isNonnull()) {
            has_kf = true;
            dpt_kf=(RefKF->ptError()*RefKF->ptError());
            nhit_kf=RefKF->hitPattern().trackerLayersWithMeasurement();
            momentum_kf.SetPx(RefKF->px());
            momentum_kf.SetPy(RefKF->py());
            momentum_kf.SetPz(RefKF->pz());
            float ENE=sqrt(RefKF->p()*RefKF->p()+m_el*m_el);
            if( DebugIDCandidates ) 
              cout << "SetCandidates:: KFTrackRef: Ene " << ENE << endl;
            
            momentum_kf.SetE(ENE);
          }
          else {
            if( DebugIDCandidates ) 
              cout <<  "SetCandidates:: !!!! NULL KF Track Ref " << endl;
          }
        } 
        else {
          if( DebugIDCandidates ) 
            cout <<  "SetCandidates:: !!!! NotAvailable KF Track Ref " << endl;
        }
      } 

      vector<double> ps1Ene(0);
      vector<double> ps2Ene(0);
      // Important is the PS clusters are not saved before the ecal one, these
      // energy are not correctly assigned 
      // For the moment I get only the closest PS clusters: this has to be changed
      if  (typeassoele == reco::PFBlockElement::PS1) {  
        PFClusterRef  psref = elements[(index_assogsf[ielegsf])].clusterRef();
        ps1Ene.push_back(psref->energy());
        ps1TotEne += psref->energy();
        elementsToAdd.push_back((index_assogsf[ielegsf]));
      }
      if  (typeassoele == reco::PFBlockElement::PS2) {  
        PFClusterRef  psref = elements[(index_assogsf[ielegsf])].clusterRef();
        ps2Ene.push_back(psref->energy());
        ps2TotEne += psref->energy();
        elementsToAdd.push_back((index_assogsf[ielegsf]));
      }

      if  (typeassoele == reco::PFBlockElement::ECAL) {
        elementsToAdd.push_back((index_assogsf[ielegsf]));

        GSFClus_index = index_assogsf[ielegsf];

        const reco::PFBlockElementCluster * clust =  
          dynamic_cast<const reco::PFBlockElementCluster*>((&elements[(index_assogsf[ielegsf])])); 
        
        eecal++;
        
        reco::PFCluster cl=*clust->clusterRef();
        
        float EE=pfcalib_.energyEm(cl,ps1Ene,ps2Ene);
        
        float ceta=clust->clusterRef()->position().eta();
        float cphi=clust->clusterRef()->position().phi();
        
        float mphi=-2.97025;
        if (ceta<0) mphi+=0.00638;
        for (int ip=1; ip<19; ip++){
          float df= cphi - (mphi+(ip*6.283185/18));
          if (fabs(df)<0.01) goodphi=false;
        }
        float dE=pfresol_.getEnergyResolutionEm(EE,clust->clusterRef()->position().eta());
        if( DebugIDCandidates ) 
          cout << "SetCandidates:: EcalCluster: Ene " << EE << " dE " <<  dE <<endl;
        
        Eene+=EE;
        dene+=dE*dE;
      }
      
      if  (typeassoele == reco::PFBlockElement::HCAL) {
        
        const reco::PFBlockElementCluster * clust =  
          dynamic_cast<const reco::PFBlockElementCluster*>((&elements[(index_assogsf[ielegsf])])); 
        float df=fabs(momentum_gsf.phi()-clust->clusterRef()->position().phi());
        if (df<0.1) {
          // Not used for the moments
          // elementsToAdd.push_back((index_assogsf[ielegsf])); 
          //Eene+=clust->clusterRef()->energy();
          Hene+=clust->clusterRef()->energy();
          hcal++;
        }
      }

      // Important: Add energy from the brems
      if  (typeassoele == reco::PFBlockElement::BREM) {
        for (uint ibrem = 0; ibrem < associatedToBrems_.size(); ibrem++){
          vector<unsigned> index_assobrem = (associatedToBrems_[ibrem]).second;
          unsigned index_brem = (associatedToBrems_[ibrem]).first;
          if (index_brem == index_assogsf[ielegsf]) {   // Note: important: check that 
            // brem emission is from the considered gsf track

            // if the brem info are needed:
            //const reco::PFBlockElementBrem * BremEl  =  
            //  dynamic_cast<const reco::PFBlockElementBrem*>((&elements[index_brem]));
            //unsigned BremTrajP = BremEl->indTrajPoint();
            //unsigned TrajPos = BremTrajP-2;
            // if (TrajPos <= 3) do something:  This could be important because
            // in case of early brem the pin is badly estimated    
            elementsToAdd.push_back(index_brem);
            for (unsigned ielebrem=0;ielebrem<index_assobrem.size();ielebrem++) {
              vector<double> ps1EneFromBrem(0);
              vector<double> ps2EneFromBrem(0);
              // Important is the PS clusters are not saved before the ecal one, these
              // energy are not correctly assigned 
              // For the moment I get only the closest PS clusters: this has to be changed
              if (elements[(index_assobrem[ielebrem])].type() == reco::PFBlockElement::PS1) {
                PFClusterRef  psref = elements[(index_assobrem[ielebrem])].clusterRef();
                ps1EneFromBrem.push_back(psref->energy());
                ps1TotEne += psref->energy();
                elementsToAdd.push_back(index_assobrem[ielebrem]);
              }
              if (elements[(index_assobrem[ielebrem])].type() == reco::PFBlockElement::PS2) {
                PFClusterRef  psref = elements[(index_assobrem[ielebrem])].clusterRef();
                ps2EneFromBrem.push_back(psref->energy());
                ps2TotEne += psref->energy();
                elementsToAdd.push_back(index_assobrem[ielebrem]);
              }        

              if (elements[(index_assobrem[ielebrem])].type() == reco::PFBlockElement::ECAL) {
                if( index_assobrem[ielebrem] !=  GSFClus_index) {
                  elementsToAdd.push_back(index_assobrem[ielebrem]);
                  reco::PFClusterRef clusterRef = elements[(index_assobrem[ielebrem])].clusterRef();
                  
                  float EE = pfcalib_.energyEm(*clusterRef,ps1EneFromBrem,ps2EneFromBrem);
                  float ceta = clusterRef->position().eta();
                  // float cphi = clusterRef->position().phi();
                  float dE=pfresol_.getEnergyResolutionEm(EE,ceta);
                  if( DebugIDCandidates ) cout << "SetCandidates:: BremCluster: Ene " << EE << " dE " <<  dE <<endl;      
                  Eene+=EE;
                  dene+=dE*dE;
                }
              } 
              if (elements[(index_assobrem[ielebrem])].type() == reco::PFBlockElement::HCAL) {
                reco::PFClusterRef clusterRef = elements[(index_assobrem[ielebrem])].clusterRef();
                float df=fabs(momentum_gsf.phi()-clusterRef->position().phi());
                if (df<0.1) {
                  // Not used for the moments
                  // elementsToAdd.push_back(index_assobrem[ielebrem]); 
                  //Eene+= clusterRef->energy();
                  Hene+= clusterRef->energy();
                  hcal++;
                }
              }
            }
          }
        }
      }
    } // End Loop On element associated to the GSF tracks
    if (has_gsf) {
      if ((nhit_gsf<8) && (has_kf)){
        
        // Use Hene if some condition.... 
        
        
        de_gs=1-momentum_gsf.E()/Eene;
        de_me=1-momentum_mean.E()/Eene;
        de_kf=1-momentum_kf.E()/Eene;
        
        
        momentum=momentum_kf;
        float Fe=Eene;
        
        float scale= Fe/momentum.E();
        
        newmomentum.SetPxPyPzE(scale*momentum.Px(),
                               scale*momentum.Py(),
                               scale*momentum.Pz(),Fe);
        //      isvalid_=true;
        
      } 
      if ((nhit_gsf>7) || (has_kf==false)){
        
        de_gs=1-momentum_gsf.E()/Eene;
        de_me=1-momentum_mean.E()/Eene;
        de_kf=1-momentum_kf.E()/Eene;
        
        momentum=momentum_gsf;
        dpt=1/(dpt_gsf*dpt_gsf);
        
        dene= 1/dene;
        
        float Fe=((dene*Eene) +(dpt*momentum.E()))/(dene+dpt);
        
        if ((de_gs>0.05)&&(de_kf>0.05)){
          Fe=Eene;
        }
        if ((de_gs<-0.1)&&(de_me<-0.1) &&(de_kf<0.)){
          Fe=momentum.E();
        }
        float scale= Fe/momentum.E();
        
        newmomentum.SetPxPyPzE(scale*momentum.Px(),
                               scale*momentum.Py(),
                               scale*momentum.Pz(),Fe);
        
      }
      if (newmomentum.pt()>0.5){
        
        // the pf candidate are created: we need to set something more? 
        // IMPORTANT -> We need the gsftrackRef, not only the TrackRef??

        if( DebugIDCandidates )
          cout << "SetCandidates:: I am before doing candidate " <<endl;
        
        reco::PFCandidate::ParticleType particleType 
          = reco::PFCandidate::e;
        reco::PFCandidate temp_Candidate;
        temp_Candidate = PFCandidate(charge,newmomentum,particleType);
        temp_Candidate.set_mva_e_pi(BDToutput_[igsf]);
        temp_Candidate.setEcalEnergy(Eene);
        temp_Candidate.setHcalEnergy(0.);  // Use HCAL energy? 
        temp_Candidate.setPs1Energy(ps1TotEne);
        temp_Candidate.setPs2Energy(ps2TotEne);
        temp_Candidate.setTrackRef(RefKF);   
        // This reference could be NULL it is needed a protection? 
        temp_Candidate.setGsfTrackRef(RefGSF);
        
        if( DebugIDCandidates ) 
          cout << "SetCandidates:: I am after doing candidate " <<endl;
        
        for (uint elad=0; elad<elementsToAdd.size();elad++){
          temp_Candidate.addElementInBlock(blockRef,elementsToAdd[elad]);
        }
        
        elCandidate_.push_back(temp_Candidate);

      }
      else {
        BDToutput_[igsf] = -1.;   // if the momentum is < 0.5 ID = false, but not sure
        // it could be misleading. 
        if( DebugIDCandidates ) 
          cout << "SetCandidates:: No Candidate Produced because of Pt cut: 0.5 " <<endl;
      }
    } 
    else {
      BDToutput_[igsf] = -1.;  // if gsf ref does not exist
      if( DebugIDCandidates ) 
        cout << "SetCandidates:: No Candidate Produced because of No GSF Track Ref " <<endl;
    }
  } // End Loop On GSF tracks
  return;
}
// This function get the associatedToGSF and associatedToBrems maps and for each pfelectron candidate 
// the element used are set active == false. 
// note: the HCAL cluster are not used for the moments and also are not locked. 
void PFElectronAlgo::SetActive(const reco::PFBlockRef&  blockRef,
                               const AssMap& associatedToGSF_,
                               const AssMap& associatedToBrems_,
                               std::vector<bool>& active){
  const reco::PFBlock& block = *blockRef;
  PFBlock::LinkData linkData =  block.linkData();     
  const edm::OwnVector< reco::PFBlockElement >&  elements = block.elements();
  for (uint igsf = 0; igsf<associatedToGSF_.size();igsf++) {
    if(BDToutput_[igsf] < mvaEleCut_) continue;  // Important enter only if good cut: to set configurable
    uint index_gsf =  associatedToGSF_[igsf].first;
    active[index_gsf] = false;  // the gsf
    vector<uint> index_assogsf =  associatedToGSF_[igsf].second;
    for  (unsigned ielegsf=0;ielegsf<index_assogsf.size();ielegsf++) {
      PFBlockElement::Type typeassoele = elements[(index_assogsf[ielegsf])].type();
      if (typeassoele == reco::PFBlockElement::HCAL) continue;  // the HCAL elements are not disattivated 
      active[(index_assogsf[ielegsf])] = false;  // the elements associated to the gsf
      if (typeassoele == reco::PFBlockElement::BREM) {
        for (uint ibrem = 0; ibrem < associatedToBrems_.size(); ibrem++){
          vector<unsigned> index_assobrem = (associatedToBrems_[ibrem]).second;
          unsigned index_brem = (associatedToBrems_[ibrem]).first;       
          if (index_brem == index_assogsf[ielegsf]) {   // Note: important: check that 
            // brem emission is from the considered gsf track
            for (unsigned ielebrem=0;ielebrem<index_assobrem.size();ielebrem++) {
              if (elements[(index_assobrem[ielebrem])].type() == reco::PFBlockElement::HCAL) continue;
              active[(index_assobrem[ielebrem])] = false;  // the elements associated to the brems
            }
          }
        }  // End loop on elements from brem      
      }
    } // End loop on elements from gsf track
  }  // End loop on gsf track  
  return;
}

---