/data/refman/pasoursint/CMSSW_6_1_2_SLHC4_patch1/src/RecoParticleFlow/PFProducer/src/PFAlgo.cc

Go to the documentation of this file.

#include "FWCore/MessageLogger/interface/MessageLogger.h"

#include "RecoParticleFlow/PFProducer/interface/PFAlgo.h"
#include "RecoParticleFlow/PFProducer/interface/PFMuonAlgo.h"  //PFMuons
#include "RecoParticleFlow/PFProducer/interface/PFElectronAlgo.h"  
#include "RecoParticleFlow/PFProducer/interface/PFPhotonAlgo.h"    
#include "RecoParticleFlow/PFProducer/interface/PFElectronExtraEqual.h"

#include "RecoParticleFlow/PFClusterTools/interface/PFEnergyCalibration.h"
#include "RecoParticleFlow/PFClusterTools/interface/PFEnergyCalibrationHF.h"
#include "RecoParticleFlow/PFClusterTools/interface/PFSCEnergyCalibration.h"

#include "DataFormats/ParticleFlowReco/interface/PFRecHit.h"
#include "DataFormats/ParticleFlowReco/interface/PFBlock.h"
#include "DataFormats/ParticleFlowReco/interface/PFBlockElementTrack.h"
#include "DataFormats/ParticleFlowReco/interface/PFBlockElementCluster.h"
#include "DataFormats/ParticleFlowReco/interface/PFRecTrack.h"
#include "DataFormats/ParticleFlowReco/interface/PFCluster.h"
#include "DataFormats/ParticleFlowReco/interface/PFLayer.h"

#include "DataFormats/VertexReco/interface/Vertex.h"
#include "DataFormats/TrackReco/interface/Track.h"
#include "DataFormats/MuonReco/interface/Muon.h"

// #include "FWCore/Framework/interface/OrphanHandle.h"
#include "DataFormats/Common/interface/OrphanHandle.h"
// #include "DataFormats/Common/interface/ProductID.h"
#include "DataFormats/Provenance/interface/ProductID.h"
#include "DataFormats/Math/interface/LorentzVector.h"


#include "DataFormats/ParticleFlowCandidate/interface/PFCandidate.h"
#include "DataFormats/ParticleFlowCandidate/interface/PFCandidateFwd.h"

#include "DataFormats/ParticleFlowReco/interface/PFDisplacedVertex.h"
#include "DataFormats/ParticleFlowReco/interface/PFDisplacedVertexFwd.h"


#include "Math/PxPyPzM4D.h"
#include "Math/LorentzVector.h"
#include "Math/DisplacementVector3D.h"
#include "Math/SMatrix.h"
#include "TDecompChol.h"

#include "boost/graph/adjacency_matrix.hpp" 
#include "boost/graph/graph_utility.hpp" 


using namespace std;
using namespace reco;
using namespace boost;

typedef std::list< reco::PFBlockRef >::iterator IBR;



PFAlgo::PFAlgo()
  : pfCandidates_( new PFCandidateCollection),
    nSigmaECAL_(0), 
    nSigmaHCAL_(1),
    algo_(1),
    debug_(false),
    pfele_(0),
    pfpho_(0),
    useVertices_(false)
{}

PFAlgo::~PFAlgo() {
  if (usePFElectrons_) delete pfele_;
  if (usePFPhotons_)     delete pfpho_;
}


void
PFAlgo::setParameters(double nSigmaECAL,
                      double nSigmaHCAL, 
                      const boost::shared_ptr<PFEnergyCalibration>& calibration,
                      const boost::shared_ptr<PFEnergyCalibrationHF>&  thepfEnergyCalibrationHF) {

  nSigmaECAL_ = nSigmaECAL;
  nSigmaHCAL_ = nSigmaHCAL;

  calibration_ = calibration;
  thepfEnergyCalibrationHF_ = thepfEnergyCalibrationHF;

}

//PFElectrons: a new method added to set the parameters for electron reconstruction. 
void 
PFAlgo::setPFEleParameters(double mvaEleCut,
                           string mvaWeightFileEleID,
                           bool usePFElectrons,
                           const boost::shared_ptr<PFSCEnergyCalibration>& thePFSCEnergyCalibration,
                           const boost::shared_ptr<PFEnergyCalibration>& thePFEnergyCalibration,
                           double sumEtEcalIsoForEgammaSC_barrel,
                           double sumEtEcalIsoForEgammaSC_endcap,
                           double coneEcalIsoForEgammaSC,
                           double sumPtTrackIsoForEgammaSC_barrel,
                           double sumPtTrackIsoForEgammaSC_endcap,
                           unsigned int nTrackIsoForEgammaSC,
                           double coneTrackIsoForEgammaSC,
                           bool applyCrackCorrections,
                           bool usePFSCEleCalib,
                           bool useEGElectrons,
                           bool useEGammaSupercluster) {
  
  mvaEleCut_ = mvaEleCut;
  usePFElectrons_ = usePFElectrons;
  applyCrackCorrectionsElectrons_ = applyCrackCorrections;  
  usePFSCEleCalib_ = usePFSCEleCalib;
  thePFSCEnergyCalibration_ = thePFSCEnergyCalibration;
  useEGElectrons_ = useEGElectrons;
  useEGammaSupercluster_ = useEGammaSupercluster;
  sumEtEcalIsoForEgammaSC_barrel_ = sumEtEcalIsoForEgammaSC_barrel;
  sumEtEcalIsoForEgammaSC_endcap_ = sumEtEcalIsoForEgammaSC_endcap;
  coneEcalIsoForEgammaSC_ = coneEcalIsoForEgammaSC;
  sumPtTrackIsoForEgammaSC_barrel_ = sumPtTrackIsoForEgammaSC_barrel;
  sumPtTrackIsoForEgammaSC_endcap_ = sumPtTrackIsoForEgammaSC_endcap;
  coneTrackIsoForEgammaSC_ = coneTrackIsoForEgammaSC;
  nTrackIsoForEgammaSC_ = nTrackIsoForEgammaSC;


  if(!usePFElectrons_) return;
  mvaWeightFileEleID_ = mvaWeightFileEleID;
  FILE * fileEleID = fopen(mvaWeightFileEleID_.c_str(), "r");
  if (fileEleID) {
    fclose(fileEleID);
  }
  else {
    string err = "PFAlgo: cannot open weight file '";
    err += mvaWeightFileEleID;
    err += "'";
    throw invalid_argument( err );
  }
  pfele_= new PFElectronAlgo(mvaEleCut_,mvaWeightFileEleID_,
                             thePFSCEnergyCalibration_,
                             thePFEnergyCalibration,
                             applyCrackCorrectionsElectrons_,
                             usePFSCEleCalib_,
                             useEGElectrons_,
                             useEGammaSupercluster_,
                             sumEtEcalIsoForEgammaSC_barrel_,
                             sumEtEcalIsoForEgammaSC_endcap_,
                             coneEcalIsoForEgammaSC_,
                             sumPtTrackIsoForEgammaSC_barrel_,
                             sumPtTrackIsoForEgammaSC_endcap_,
                             nTrackIsoForEgammaSC_,
                             coneTrackIsoForEgammaSC_);
}

void 
PFAlgo::setPFPhotonParameters(bool usePFPhotons,  
                              std::string mvaWeightFileConvID, 
                              double mvaConvCut,
                              bool useReg,
                              std::string X0_Map,
                              const boost::shared_ptr<PFEnergyCalibration>& thePFEnergyCalibration,
                              double sumPtTrackIsoForPhoton,
                              double sumPtTrackIsoSlopeForPhoton)
 {

  usePFPhotons_ = usePFPhotons;

  //for MVA pass PV if there is one in the collection otherwise pass a dummy    
  reco::Vertex dummy;  
  if(useVertices_)  
    {  
      dummy = primaryVertex_;  
    }  
  else { // create a dummy PV  
    reco::Vertex::Error e;  
    e(0, 0) = 0.0015 * 0.0015;  
    e(1, 1) = 0.0015 * 0.0015;  
    e(2, 2) = 15. * 15.;  
    reco::Vertex::Point p(0, 0, 0);  
    dummy = reco::Vertex(p, e, 0, 0, 0);  
  }  
  // pv=&dummy;  
  if(! usePFPhotons_) return;  
  FILE * filePhotonConvID = fopen(mvaWeightFileConvID.c_str(), "r");  
  if (filePhotonConvID) {  
    fclose(filePhotonConvID);  
  }  
  else {  
    string err = "PFAlgo: cannot open weight file '";  
    err += mvaWeightFileConvID;  
    err += "'";  
    throw invalid_argument( err );  
  }  
  const reco::Vertex* pv=&dummy;  
  pfpho_ = new PFPhotonAlgo(mvaWeightFileConvID, 
                            mvaConvCut, 
                            useReg,
                            X0_Map,  
                            *pv,
                            thePFEnergyCalibration,
                            sumPtTrackIsoForPhoton,
                            sumPtTrackIsoSlopeForPhoton
                            );
  return;
}
/*
void PFAlgo::setPFPhotonRegWeights(
                  const GBRForest *LCorrForest,
                  const GBRForest *GCorrForest,
                  const GBRForest *ResForest
                  ) {                                                           
  if(usePFPhotons_) 
    pfpho_->setGBRForest(LCorrForest, GCorrForest, ResForest);
} 
*/
void PFAlgo::setPFPhotonRegWeights(
                                   const GBRForest *LCorrForestEB,
                                   const GBRForest *LCorrForestEE,
                                   const GBRForest *GCorrForestBarrel,
                                   const GBRForest *GCorrForestEndcapHr9,
                                   const GBRForest *GCorrForestEndcapLr9,                                          const GBRForest *PFEcalResolution
                                   ){
  
  pfpho_->setGBRForest(LCorrForestEB,LCorrForestEE,
                       GCorrForestBarrel, GCorrForestEndcapHr9, 
                       GCorrForestEndcapLr9, PFEcalResolution);
}
void 
PFAlgo::setPFMuonAndFakeParameters(std::vector<double> muonHCAL,
                                   std::vector<double> muonECAL,
                                   std::vector<double> muonHO,
                                   double nSigmaTRACK,
                                   double ptError,
                                   std::vector<double> factors45,
                                   bool usePFMuonMomAssign,
                                   bool useBestMuonTrack) 
{
  muonHCAL_ = muonHCAL;
  muonECAL_ = muonECAL;
  muonHO_ = muonHO;
  nSigmaTRACK_ = nSigmaTRACK;
  ptError_ = ptError;
  factors45_ = factors45;
  usePFMuonMomAssign_ = usePFMuonMomAssign;
  useBestMuonTrack_ = useBestMuonTrack;
}
  
void 
PFAlgo::setPostHFCleaningParameters(bool postHFCleaning,
                                    double minHFCleaningPt,
                                    double minSignificance,
                                    double maxSignificance,
                                    double minSignificanceReduction,
                                    double maxDeltaPhiPt,
                                    double minDeltaMet) { 
  postHFCleaning_ = postHFCleaning;
  minHFCleaningPt_ = minHFCleaningPt;
  minSignificance_ = minSignificance;
  maxSignificance_ = maxSignificance;
  minSignificanceReduction_= minSignificanceReduction;
  maxDeltaPhiPt_ = maxDeltaPhiPt;
  minDeltaMet_ = minDeltaMet;
}

void 
PFAlgo::setDisplacedVerticesParameters(bool rejectTracks_Bad,
                                       bool rejectTracks_Step45,
                                       bool usePFNuclearInteractions,
                                       bool usePFConversions,
                                       bool usePFDecays, 
                                       double dptRel_DispVtx){

  rejectTracks_Bad_ = rejectTracks_Bad;
  rejectTracks_Step45_ = rejectTracks_Step45;
  usePFNuclearInteractions_ = usePFNuclearInteractions;
  usePFConversions_ = usePFConversions;
  usePFDecays_ = usePFDecays;
  dptRel_DispVtx_ = dptRel_DispVtx;

}


void
PFAlgo::setPFVertexParameters(bool useVertex,
                              const reco::VertexCollection& primaryVertices) {
  useVertices_ = useVertex;
  //Now find the primary vertex!
  bool primaryVertexFound = false;
  int nVtx=primaryVertices.size();
  if(usePFPhotons_){
    pfpho_->setnPU(nVtx);
  }
  for (unsigned short i=0 ;i<primaryVertices.size();++i)
    {
      if(primaryVertices[i].isValid()&&(!primaryVertices[i].isFake()))
        {
          primaryVertex_ = primaryVertices[i];
          primaryVertexFound = true;
          break;
        }
    }
  //Use vertices if the user wants to but only if it exists a good vertex 
  useVertices_ = useVertex && primaryVertexFound; 
  if(usePFPhotons_) {
    if (useVertices_ ){
      pfpho_->setPhotonPrimaryVtx(primaryVertex_ );
    }
    else{
      reco::Vertex::Error e;  
      e(0, 0) = 0.0015 * 0.0015;  
      e(1, 1) = 0.0015 * 0.0015;  
      e(2, 2) = 15. * 15.;  
      reco::Vertex::Point p(0, 0, 0);  
      reco::Vertex dummy = reco::Vertex(p, e, 0, 0, 0); 
      //      std::cout << " PFPho " << pfpho_ << std::endl;
      pfpho_->setPhotonPrimaryVtx(dummy);
    }
  }
}


void PFAlgo::reconstructParticles( const reco::PFBlockHandle& blockHandle ) {

  blockHandle_ = blockHandle;
  reconstructParticles( *blockHandle_ ); 
}



void PFAlgo::reconstructParticles( const reco::PFBlockCollection& blocks ) {

  // reset output collection
  if(pfCandidates_.get() )
    pfCandidates_->clear();
  else 
    pfCandidates_.reset( new reco::PFCandidateCollection );

  if(pfElectronCandidates_.get() )
    pfElectronCandidates_->clear();
  else
    pfElectronCandidates_.reset( new reco::PFCandidateCollection);

  // Clearing pfPhotonCandidates
  if( pfPhotonCandidates_.get() )
    pfPhotonCandidates_->clear();
  else
    pfPhotonCandidates_.reset( new reco::PFCandidateCollection);

  if(pfCleanedCandidates_.get() ) 
    pfCleanedCandidates_->clear();
  else
    pfCleanedCandidates_.reset( new reco::PFCandidateCollection );
  
  if(pfCosmicsMuonCleanedCandidates_.get() )
    pfCosmicsMuonCleanedCandidates_->clear();
  else 
    pfCosmicsMuonCleanedCandidates_.reset( new reco::PFCandidateCollection );

  if(pfCleanedTrackerAndGlobalMuonCandidates_.get() )
    pfCleanedTrackerAndGlobalMuonCandidates_->clear();
  else 
    pfCleanedTrackerAndGlobalMuonCandidates_.reset( new reco::PFCandidateCollection );

  if(pfFakeMuonCleanedCandidates_.get() )
    pfFakeMuonCleanedCandidates_->clear();
  else 
    pfFakeMuonCleanedCandidates_.reset( new reco::PFCandidateCollection );

  if(pfPunchThroughMuonCleanedCandidates_.get() )
    pfPunchThroughMuonCleanedCandidates_->clear();
  else 
    pfPunchThroughMuonCleanedCandidates_.reset( new reco::PFCandidateCollection );

  if(pfPunchThroughHadronCleanedCandidates_.get() )
    pfPunchThroughHadronCleanedCandidates_->clear();
  else 
    pfPunchThroughHadronCleanedCandidates_.reset( new reco::PFCandidateCollection );

  if(pfAddedMuonCandidates_.get() )
    pfAddedMuonCandidates_->clear();
  else 
    pfAddedMuonCandidates_.reset( new reco::PFCandidateCollection );

  // not a auto_ptr; shout not be deleted after transfer
  pfElectronExtra_.clear();
  pfPhotonExtra_.clear();
  
  if( debug_ ) {
    cout<<"*********************************************************"<<endl;
    cout<<"*****           Particle flow algorithm             *****"<<endl;
    cout<<"*********************************************************"<<endl;
  }

  // sort elements in three lists:
  std::list< reco::PFBlockRef > hcalBlockRefs;
  std::list< reco::PFBlockRef > ecalBlockRefs;
  std::list< reco::PFBlockRef > hoBlockRefs;
  std::list< reco::PFBlockRef > otherBlockRefs;
  
  for( unsigned i=0; i<blocks.size(); ++i ) {
    // reco::PFBlockRef blockref( blockh,i );
    reco::PFBlockRef blockref = createBlockRef( blocks, i);
    
    const reco::PFBlock& block = *blockref;
    const edm::OwnVector< reco::PFBlockElement >& 
      elements = block.elements();
   
    bool singleEcalOrHcal = false;
    if( elements.size() == 1 ){
      if( elements[0].type() == reco::PFBlockElement::ECAL ){
        ecalBlockRefs.push_back( blockref );
        singleEcalOrHcal = true;
      }
      if( elements[0].type() == reco::PFBlockElement::HCAL ){
        hcalBlockRefs.push_back( blockref );
        singleEcalOrHcal = true;
      }
      if( elements[0].type() == reco::PFBlockElement::HO ){
        // Single HO elements are likely to be noise. Not considered for now.
        hoBlockRefs.push_back( blockref );
        singleEcalOrHcal = true;
      }
    }
    
    if(!singleEcalOrHcal) {
      otherBlockRefs.push_back( blockref );
    }
  }//loop blocks
  
  if( debug_ ){
    cout<<"# Ecal blocks: "<<ecalBlockRefs.size()
        <<", # Hcal blocks: "<<hcalBlockRefs.size()
        <<", # HO blocks: "<<hoBlockRefs.size()
        <<", # Other blocks: "<<otherBlockRefs.size()<<endl;
  }


  // loop on blocks that are not single ecal, 
  // and not single hcal.

  unsigned nblcks = 0;
  for( IBR io = otherBlockRefs.begin(); io!=otherBlockRefs.end(); ++io) {
    if ( debug_ ) std::cout << "Block number " << nblcks++ << std::endl; 
    processBlock( *io, hcalBlockRefs, ecalBlockRefs );
  }

  std::list< reco::PFBlockRef > empty;

  unsigned hblcks = 0;
  // process remaining single hcal blocks
  for( IBR ih = hcalBlockRefs.begin(); ih!=hcalBlockRefs.end(); ++ih) {
    if ( debug_ ) std::cout << "HCAL block number " << hblcks++ << std::endl;
    processBlock( *ih, empty, empty );
  }

  unsigned eblcks = 0;
  // process remaining single ecal blocks
  for( IBR ie = ecalBlockRefs.begin(); ie!=ecalBlockRefs.end(); ++ie) {
    if ( debug_ ) std::cout << "ECAL block number " << eblcks++ << std::endl;
    processBlock( *ie, empty, empty );
  }

  // Post HF Cleaning
  postCleaning();
}


void PFAlgo::processBlock( const reco::PFBlockRef& blockref,
                           std::list<reco::PFBlockRef>& hcalBlockRefs, 
                           std::list<reco::PFBlockRef>& ecalBlockRefs ) { 
  
  // debug_ = false;
  assert(!blockref.isNull() );
  const reco::PFBlock& block = *blockref;

  typedef std::multimap<double, unsigned>::iterator IE;
  typedef std::multimap<double, std::pair<unsigned,math::XYZVector> >::iterator IS;
  typedef std::multimap<double, std::pair<unsigned,bool> >::iterator IT;
  typedef std::multimap< unsigned, std::pair<double, unsigned> >::iterator II;

  if(debug_) {
    cout<<"#########################################################"<<endl;
    cout<<"#####           Process Block:                      #####"<<endl;
    cout<<"#########################################################"<<endl;
    cout<<block<<endl;
  }


  const edm::OwnVector< reco::PFBlockElement >& elements = block.elements();
  // make a copy of the link data, which will be edited.
  PFBlock::LinkData linkData =  block.linkData();
  
  // keep track of the elements which are still active.
  vector<bool>   active( elements.size(), true );


  // //PFElectrons:
  // usePFElectrons_ external configurable parameter to set the usage of pf electron
  std::vector<reco::PFCandidate> tempElectronCandidates;
  tempElectronCandidates.clear();
  if (usePFElectrons_) {
    if (pfele_->isElectronValidCandidate(blockref,active, primaryVertex_ )){
      // if there is at least a valid candidate is get the vector of pfcandidates
      const std::vector<reco::PFCandidate> PFElectCandidates_(pfele_->getElectronCandidates());
      for ( std::vector<reco::PFCandidate>::const_iterator ec=PFElectCandidates_.begin();   ec != PFElectCandidates_.end(); ++ec )tempElectronCandidates.push_back(*ec);
      
      // (***) We're filling the ElectronCandidates into the PFCandiate collection
      // ..... Once we let PFPhotonAlgo over-write electron-decision, we need to move this to
      // ..... after the PhotonAlgo has run (Fabian)
    }    
    // The vector active is automatically changed (it is passed by ref) in PFElectronAlgo
    // for all the electron candidate      
    pfElectronCandidates_->insert(pfElectronCandidates_->end(), 
                                  pfele_->getAllElectronCandidates().begin(), 
                                  pfele_->getAllElectronCandidates().end()); 

    pfElectronExtra_.insert(pfElectronExtra_.end(), 
                            pfele_->getElectronExtra().begin(),
                            pfele_->getElectronExtra().end());
    
  }
  if( /* --- */ usePFPhotons_ /* --- */ ) {    
    
    if(debug_)
      cout<<endl<<"--------------- entering PFPhotonAlgo ----------------"<<endl;
    vector<PFCandidatePhotonExtra> pfPhotonExtraCand;
    if ( pfpho_->isPhotonValidCandidate(blockref,               // passing the reference to the PFBlock
                                        active,                 // std::vector<bool> containing information about acitivity
                                        pfPhotonCandidates_,    // pointer to candidate vector, to be filled by the routine
                                        pfPhotonExtraCand,      // candidate extra vector, to be filled by the routine   
                                        tempElectronCandidates
                                        //pfElectronCandidates_   // pointer to some auziliary UNTOUCHED FOR NOW
                                        ) ) {
      if(debug_)
        std::cout<< " In this PFBlock we found "<<pfPhotonCandidates_->size()<<" Photon Candidates."<<std::endl;
      
      // CAUTION: In case we want to allow the PhotonAlgo to 'over-write' what the ElectronAlgo did above
      // ........ we should NOT fill the PFCandidate-vector with the electrons above (***)
      
      // Here we need to add all the photon cands to the pfCandidate list
      unsigned int extracand =0;
      PFCandidateCollection::const_iterator cand = pfPhotonCandidates_->begin();      
      for( ; cand != pfPhotonCandidates_->end(); ++cand, ++extracand) {
        pfCandidates_->push_back(*cand);
        pfPhotonExtra_.push_back(pfPhotonExtraCand[extracand]);
      }
      
    } // end of 'if' in case photons are found    
    pfPhotonExtraCand.clear();
    pfPhotonCandidates_->clear();
  } // end of Photon algo
  
  //Lock extra conversion tracks not used by Photon Algo
  if (usePFConversions_  )
    {
      for(unsigned iEle=0; iEle<elements.size(); iEle++) {
        PFBlockElement::Type type = elements[iEle].type();
        if(type==PFBlockElement::TRACK)
          {
            if(elements[iEle].trackRef()->algo() == 12)
              active[iEle]=false;       
            if(elements[iEle].trackRef()->quality(reco::TrackBase::highPurity))continue;
            const reco::PFBlockElementTrack * trackRef = dynamic_cast<const reco::PFBlockElementTrack*>((&elements[iEle]));
            if(!(trackRef->trackType(reco::PFBlockElement::T_FROM_GAMMACONV)))continue;
            if(elements[iEle].convRef().isNonnull())active[iEle]=false;
          }
      }
    }
  for ( std::vector<reco::PFCandidate>::const_iterator ec=tempElectronCandidates.begin();   ec != tempElectronCandidates.end(); ++ec ){
    pfCandidates_->push_back(*ec);  
  } 
  tempElectronCandidates.clear();
  
  
  if(debug_) 
    cout<<endl<<"--------------- loop 1 ------------------"<<endl;

  //COLINFEB16
  // In loop 1, we loop on all elements. 

  // The primary goal is to deal with tracks that are:
  // - not associated to an HCAL cluster
  // - not identified as an electron. 
  // Those tracks should be predominantly relatively low energy charged 
  // hadrons which are not detected in the ECAL. 

  // The secondary goal is to prepare for the next loops
  // - The ecal and hcal elements are sorted in separate vectors
  // which will be used as a base for the corresponding loops.
  // - For tracks which are connected to more than one HCAL cluster, 
  // the links between the track and the cluster are cut for all clusters 
  // but the closest one. 
  // - HF only blocks ( HFEM, HFHAD, HFEM+HFAD) are identified 

  // obsolete comments? 
  // loop1: 
  // - sort ecal and hcal elements in separate vectors
  // - for tracks:
  //       - lock closest ecal cluster 
  //       - cut link to farthest hcal cluster, if more than 1.

  // vectors to store indices to ho, hcal and ecal elements 
  vector<unsigned> hcalIs;
  vector<unsigned> hoIs;
  vector<unsigned> ecalIs;
  vector<unsigned> trackIs;
  vector<unsigned> ps1Is;
  vector<unsigned> ps2Is;

  vector<unsigned> hfEmIs;
  vector<unsigned> hfHadIs;
  

  for(unsigned iEle=0; iEle<elements.size(); iEle++) {
    PFBlockElement::Type type = elements[iEle].type();

    if(debug_ && type != PFBlockElement::BREM ) cout<<endl<<elements[iEle];   

    switch( type ) {
    case PFBlockElement::TRACK:
      if ( active[iEle] ) { 
        trackIs.push_back( iEle );
        if(debug_) cout<<"TRACK, stored index, continue"<<endl;
      }
      break;
    case PFBlockElement::ECAL: 
      if ( active[iEle]  ) { 
        ecalIs.push_back( iEle );
        if(debug_) cout<<"ECAL, stored index, continue"<<endl;
      }
      continue;
    case PFBlockElement::HCAL:
      if ( active[iEle] ) { 
        hcalIs.push_back( iEle );
        if(debug_) cout<<"HCAL, stored index, continue"<<endl;
      }
      continue;
    case PFBlockElement::HO:
      if (useHO_) {
        if ( active[iEle] ) { 
          hoIs.push_back( iEle );
          if(debug_) cout<<"HO, stored index, continue"<<endl;
        }
      }
      continue;
    case PFBlockElement::HFEM:
      if ( active[iEle] ) { 
        hfEmIs.push_back( iEle );
        if(debug_) cout<<"HFEM, stored index, continue"<<endl;
      }
      continue;
    case PFBlockElement::HFHAD:
      if ( active[iEle] ) { 
        hfHadIs.push_back( iEle );
        if(debug_) cout<<"HFHAD, stored index, continue"<<endl;
      }
      continue;

      //COLINFEB16 the following is not used... removing
      
      //     case PFBlockElement::PS1:
      //        if ( active[iEle] ) { 
      //         ps1Is.push_back( iEle );
      //         if(debug_) cout<<"PS1, stored index, continue"<<endl;
      //        }
      //        continue;
      //     case PFBlockElement::PS2:
      //       if ( active[iEle] ) { 
      //        ps2Is.push_back( iEle );
      //        if(debug_) cout<<"PS2, stored index, continue"<<endl;
      //       }
      //       continue;
    default:
      continue;
    }

    // we're now dealing with a track
    unsigned iTrack = iEle;
    reco::MuonRef muonRef = elements[iTrack].muonRef();    


    if (active[iTrack] && isFromSecInt(elements[iEle], "primary")){
      bool isPrimaryTrack = elements[iEle].displacedVertexRef(PFBlockElement::T_TO_DISP)->displacedVertexRef()->isTherePrimaryTracks();
      if (isPrimaryTrack) {
        if (debug_) cout << "Primary Track reconstructed alone" << endl;
        unsigned tmpi = reconstructTrack(elements[iEle]);
        (*pfCandidates_)[tmpi].addElementInBlock( blockref, iEle );
        active[iTrack] = false;
      }
    }


    if(debug_) { 
      if ( !active[iTrack] ) 
        cout << "Already used by electrons, muons, conversions" << endl;
    }
    
    // Track already used as electron, muon, conversion?
    // Added a check on the activated element
    if ( ! active[iTrack] ) continue;
  
    reco::TrackRef trackRef = elements[iTrack].trackRef();
    assert( !trackRef.isNull() ); 

    if (debug_ ) {
      cout <<"PFAlgo:processBlock "<<" "<<trackIs.size()<<" "<<ecalIs.size()<<" "<<hcalIs.size()<<" "<<hoIs.size()<<endl;
    }

    // look for associated elements of all types
    //COLINFEB16 
    // all types of links are considered. 
    // the elements are sorted by increasing distance
    std::multimap<double, unsigned> ecalElems;
    block.associatedElements( iTrack,  linkData,
                              ecalElems ,
                              reco::PFBlockElement::ECAL,
                              reco::PFBlock::LINKTEST_ALL );
    
    std::multimap<double, unsigned> hcalElems;
    block.associatedElements( iTrack,  linkData,
                              hcalElems,
                              reco::PFBlockElement::HCAL,
                              reco::PFBlock::LINKTEST_ALL );

    // When a track has no HCAL cluster linked, but another track is linked to the same
    // ECAL cluster and an HCAL cluster, link the track to the HCAL cluster for 
    // later analysis
    if ( hcalElems.empty() && !ecalElems.empty() ) {
      // debug_ = true;
      unsigned ntt = 1;
      unsigned index = ecalElems.begin()->second;
      std::multimap<double, unsigned> sortedTracks;
      block.associatedElements( index,  linkData,
                                sortedTracks,
                                reco::PFBlockElement::TRACK,
                                reco::PFBlock::LINKTEST_ALL );
      // std::cout << "The ECAL is linked to " << sortedTracks.size() << std::endl;

      // Loop over all tracks
      for(IE ie = sortedTracks.begin(); ie != sortedTracks.end(); ++ie ) {
        unsigned jTrack = ie->second;
        // std::cout << "Track " << jTrack << std::endl;
        // Track must be active
        if ( !active[jTrack] ) continue;
        //std::cout << "Active " << std::endl;
        
        // The loop is on the other tracks !
        if ( jTrack == iTrack ) continue;
        //std::cout << "A different track ! " << std::endl;
        
        // Check if the ECAL closest to this track is the current ECAL
        // Otherwise ignore this track in the neutral energy determination
        std::multimap<double, unsigned> sortedECAL;
        block.associatedElements( jTrack,  linkData,
                                  sortedECAL,
                                  reco::PFBlockElement::ECAL,
                                  reco::PFBlock::LINKTEST_ALL );
        if ( sortedECAL.begin()->second != index ) continue;
        //std::cout << "With closest ECAL identical " << std::endl;

        
        // Check if this track is also linked to an HCAL 
        std::multimap<double, unsigned> sortedHCAL;
        block.associatedElements( jTrack,  linkData,
                                  sortedHCAL,
                                  reco::PFBlockElement::HCAL,
                                  reco::PFBlock::LINKTEST_ALL );
        if ( !sortedHCAL.size() ) continue;
        //std::cout << "With an HCAL cluster " << sortedHCAL.begin()->first << std::endl;
        ntt++;
        
        // In that case establish a link with the first track 
        block.setLink( iTrack, 
                       sortedHCAL.begin()->second, 
                       sortedECAL.begin()->first, 
                       linkData,
                       PFBlock::LINKTEST_RECHIT );
        
      } // End other tracks
      
      // Redefine HCAL elements 
      block.associatedElements( iTrack,  linkData,
                                hcalElems,
                                reco::PFBlockElement::HCAL,
                                reco::PFBlock::LINKTEST_ALL );

      if ( debug_ && hcalElems.size() ) 
        std::cout << "Track linked back to HCAL due to ECAL sharing with other tracks" << std::endl;
    }
    
    //MICHELE 
    //TEMPORARY SOLUTION FOR ELECTRON REJECTION IN PFTAU
    //COLINFEB16 
    // in case particle flow electrons are not reconstructed, 
    // the mva_e_pi of the charged hadron will be set to 1 
    // if a GSF element is associated to the current TRACK element
    // This information will be used in the electron rejection for tau ID.
    std::multimap<double,unsigned> gsfElems;
    if (usePFElectrons_ == false) {
      block.associatedElements( iTrack,  linkData,
                                gsfElems ,
                                reco::PFBlockElement::GSF );
    }
    //

    if(hcalElems.empty() && debug_) {
      cout<<"no hcal element connected to track "<<iTrack<<endl;
    }

    // will now loop on associated elements ... 
    // typedef std::multimap<double, unsigned>::iterator IE;
   
    bool hcalFound = false;

    if(debug_) 
      cout<<"now looping on elements associated to the track"<<endl;
    
    // ... first on associated ECAL elements
    // Check if there is still a free ECAL for this track
    for(IE ie = ecalElems.begin(); ie != ecalElems.end(); ++ie ) {
      
      unsigned index = ie->second;
      // Sanity checks and optional printout
      PFBlockElement::Type type = elements[index].type();
      if(debug_) {
        double  dist  = ie->first;
        cout<<"\telement "<<elements[index]<<" linked with distance = "<< dist <<endl;
        if ( ! active[index] ) cout << "This ECAL is already used - skip it" << endl;      
      }
      assert( type == PFBlockElement::ECAL );

      // This ECAL is not free (taken by an electron?) - just skip it
      if ( ! active[index] ) continue;      

      // Flag ECAL clusters for which the corresponding track is not linked to an HCAL cluster

      //reco::PFClusterRef clusterRef = elements[index].clusterRef();
      //assert( !clusterRef.isNull() );
      if( !hcalElems.empty() && debug_)
        cout<<"\t\tat least one hcal element connected to the track."
            <<" Sparing Ecal cluster for the hcal loop"<<endl; 

    } //loop print ecal elements

  
    // tracks which are not linked to an HCAL 
    // are reconstructed now. 

    if( hcalElems.empty() ) {
      
      if ( debug_ ) 
        std::cout << "Now deals with tracks linked to no HCAL clusters" << std::endl; 
      // vector<unsigned> elementIndices;
      // elementIndices.push_back(iTrack); 

      // The track momentum
      double trackMomentum = elements[iTrack].trackRef()->p();      
      if ( debug_ ) 
        std::cout << elements[iTrack] << std::endl; 
      
      // Is it a "tight" muon ?
      bool thisIsAMuon = PFMuonAlgo::isMuon(elements[iTrack]);
      if ( thisIsAMuon ) trackMomentum = 0.;

      // Is it a fake track ?
      bool rejectFake = false;
      if ( !thisIsAMuon  && elements[iTrack].trackRef()->ptError() > ptError_ ) { 

        double deficit = trackMomentum; 
        double resol = neutralHadronEnergyResolution(trackMomentum,
                                                     elements[iTrack].trackRef()->eta());
        resol *= trackMomentum;

        if ( !ecalElems.empty() ) { 
          unsigned thisEcal = ecalElems.begin()->second;
          reco::PFClusterRef clusterRef = elements[thisEcal].clusterRef();
          deficit -= clusterRef->energy();
          resol = neutralHadronEnergyResolution(trackMomentum,
                                                clusterRef->positionREP().Eta());
          resol *= trackMomentum;
        }

        bool isPrimary = isFromSecInt(elements[iTrack], "primary");

        if ( deficit > nSigmaTRACK_*resol && !isPrimary) { 
          rejectFake = true;
          active[iTrack] = false;
          if ( debug_ )  
            std::cout << elements[iTrack] << std::endl
                      << "is probably a fake (1) --> lock the track" 
                      << std::endl;
        }
      }

      if ( rejectFake ) continue;

      // Create a track candidate       
      // unsigned tmpi = reconstructTrack( elements[iTrack] );
      //active[iTrack] = false;
      std::vector<unsigned> tmpi;
      std::vector<unsigned> kTrack;
      
      // Some cleaning : secondary tracks without calo's and large momentum must be fake
      double Dpt = trackRef->ptError();

      if ( rejectTracks_Step45_ && ecalElems.empty() && 
           trackMomentum > 30. && Dpt > 0.5 && 
           ( trackRef->algo() == TrackBase::iter4 || 
             trackRef->algo() == TrackBase::iter5 || 
             trackRef->algo() == TrackBase::iter6 ) ) {

        //
        double dptRel = Dpt/trackRef->pt()*100;
        bool isPrimaryOrSecondary = isFromSecInt(elements[iTrack], "all");

        if ( isPrimaryOrSecondary && dptRel < dptRel_DispVtx_) continue;
        unsigned nHits =  elements[iTrack].trackRef()->hitPattern().trackerLayersWithMeasurement();
        unsigned int NLostHit = trackRef->hitPattern().trackerLayersWithoutMeasurement();

        if ( debug_ ) 
          std::cout << "A track (algo = " << trackRef->algo() << ") with momentum " << trackMomentum 
                    << " / " << elements[iTrack].trackRef()->pt() << " +/- " << Dpt 
                    << " / " << elements[iTrack].trackRef()->eta() 
                    << " without any link to ECAL/HCAL and with " << nHits << " (" << NLostHit 
                    << ") hits (lost hits) has been cleaned" << std::endl;
        active[iTrack] = false;
        continue;
      }


      tmpi.push_back(reconstructTrack( elements[iTrack]));
      kTrack.push_back(iTrack);
      active[iTrack] = false;

      // No ECAL cluster either ... continue...
      if ( ecalElems.empty() ) { 
        (*pfCandidates_)[tmpi[0]].setEcalEnergy( 0., 0. );
        (*pfCandidates_)[tmpi[0]].setHcalEnergy( 0., 0. );
        (*pfCandidates_)[tmpi[0]].setHoEnergy( 0., 0. );
        (*pfCandidates_)[tmpi[0]].setPs1Energy( 0 );
        (*pfCandidates_)[tmpi[0]].setPs2Energy( 0 );
        (*pfCandidates_)[tmpi[0]].addElementInBlock( blockref, kTrack[0] );
        continue;
      }
          
      // Look for closest ECAL cluster
      unsigned thisEcal = ecalElems.begin()->second;
      reco::PFClusterRef clusterRef = elements[thisEcal].clusterRef();
      if ( debug_ ) std::cout << " is associated to " << elements[thisEcal] << std::endl;

      // Set ECAL energy for muons
      if ( thisIsAMuon ) { 
        (*pfCandidates_)[tmpi[0]].setEcalEnergy( clusterRef->energy(),
                                                 std::min(clusterRef->energy(), muonECAL_[0]) );
        (*pfCandidates_)[tmpi[0]].setHcalEnergy( 0., 0. );
        (*pfCandidates_)[tmpi[0]].setHoEnergy( 0., 0. );
        (*pfCandidates_)[tmpi[0]].setPs1Energy( 0 );
        (*pfCandidates_)[tmpi[0]].setPs2Energy( 0 );
        (*pfCandidates_)[tmpi[0]].addElementInBlock( blockref, kTrack[0] );
      }
      
      double slopeEcal = 1.;
      bool connectedToEcal = false;
      unsigned iEcal = 99999;
      double calibEcal = 0.;
      double calibHcal = 0.;
      double totalEcal = thisIsAMuon ? -muonECAL_[0] : 0.;

      // Consider charged particles closest to the same ECAL cluster
      std::multimap<double, unsigned> sortedTracks;
      block.associatedElements( thisEcal,  linkData,
                                sortedTracks,
                                reco::PFBlockElement::TRACK,
                                reco::PFBlock::LINKTEST_ALL );
    
      for(IE ie = sortedTracks.begin(); ie != sortedTracks.end(); ++ie ) {
        unsigned jTrack = ie->second;

        // Don't consider already used tracks
        if ( !active[jTrack] ) continue;

        // The loop is on the other tracks !
        if ( jTrack == iTrack ) continue;
        
        // Check if the ECAL cluster closest to this track is 
        // indeed the current ECAL cluster. Only tracks not 
        // linked to an HCAL cluster are considered here
        std::multimap<double, unsigned> sortedECAL;
        block.associatedElements( jTrack,  linkData,
                                  sortedECAL,
                                  reco::PFBlockElement::ECAL,
                                  reco::PFBlock::LINKTEST_ALL );
        if ( sortedECAL.begin()->second != thisEcal ) continue;

        // Check if this track is a muon
        bool thatIsAMuon = PFMuonAlgo::isMuon(elements[jTrack]);

        // Check if this track is not a fake
        bool rejectFake = false;
        reco::TrackRef trackRef = elements[jTrack].trackRef();
        if ( !thatIsAMuon && trackRef->ptError() > ptError_) { 
          double deficit = trackMomentum + trackRef->p() - clusterRef->energy();
          double resol = nSigmaTRACK_*neutralHadronEnergyResolution(trackMomentum+trackRef->p(),
                                                                    clusterRef->positionREP().Eta());
          resol *= (trackMomentum+trackRef->p());
          if ( deficit > nSigmaTRACK_*resol ) { 
            rejectFake = true;
            kTrack.push_back(jTrack);
            active[jTrack] = false;
            if ( debug_ )  
              std::cout << elements[jTrack] << std::endl
                        << "is probably a fake (2) --> lock the track" 
                        << std::endl; 
          }
        }
        if ( rejectFake ) continue;
        
        
        // Otherwise, add this track momentum to the total track momentum
        /* */
        // reco::TrackRef trackRef = elements[jTrack].trackRef();
        if ( !thatIsAMuon ) {   
          if ( debug_ ) 
            std::cout << "Track momentum increased from " << trackMomentum << " GeV "; 
          trackMomentum += trackRef->p();
          if ( debug_ ) {
            std::cout << "to " << trackMomentum << " GeV." << std::endl; 
            std::cout << "with " << elements[jTrack] << std::endl;
          }
        } else { 
          totalEcal -= muonECAL_[0];
          totalEcal = std::max(totalEcal, 0.);
        }

        // And create a charged particle candidate !
        tmpi.push_back(reconstructTrack( elements[jTrack] ));
        kTrack.push_back(jTrack);
        active[jTrack] = false;

        if ( thatIsAMuon ) { 
          (*pfCandidates_)[tmpi.back()].setEcalEnergy(clusterRef->energy(),
                                                      std::min(clusterRef->energy(),muonECAL_[0]));
          (*pfCandidates_)[tmpi.back()].setHcalEnergy( 0., 0. );
          (*pfCandidates_)[tmpi.back()].setHoEnergy( 0., 0. );
          (*pfCandidates_)[tmpi.back()].setPs1Energy( 0 );
          (*pfCandidates_)[tmpi.back()].setPs2Energy( 0 );
          (*pfCandidates_)[tmpi.back()].addElementInBlock( blockref, kTrack.back() );
        }
      }


      if ( debug_ ) std::cout << "Loop over all associated ECAL clusters" << std::endl; 
      // Loop over all ECAL linked clusters ordered by increasing distance.
      for(IE ie = ecalElems.begin(); ie != ecalElems.end(); ++ie ) {
        
        unsigned index = ie->second;
        PFBlockElement::Type type = elements[index].type();
        assert( type == PFBlockElement::ECAL );
        if ( debug_ ) std::cout << elements[index] << std::endl;
        
        // Just skip clusters already taken
        if ( debug_ && ! active[index] ) std::cout << "is not active  - ignore " << std::endl;
        if ( ! active[index] ) continue;

        // Just skip this cluster if it's closer to another track
        /* */
        block.associatedElements( index,  linkData,
                                  sortedTracks,
                                  reco::PFBlockElement::TRACK,
                                  reco::PFBlock::LINKTEST_ALL );
        bool skip = true;
        for (unsigned ic=0; ic<kTrack.size();++ic) {
          if ( sortedTracks.begin()->second == kTrack[ic] ) { 
            skip = false;
            break;
          }
        }
        if ( debug_ && skip ) std::cout << "is closer to another track - ignore " << std::endl;
        if ( skip ) continue;
        /* */
                
        // The corresponding PFCluster and energy
        reco::PFClusterRef clusterRef = elements[index].clusterRef();
        assert( !clusterRef.isNull() );

        if ( debug_ ) {
          double dist = ie->first;
          std::cout <<"Ecal cluster with raw energy = " << clusterRef->energy() 
                    <<" linked with distance = " << dist << std::endl;
        }
        /*
        double dist = ie->first;
        if ( !connectedToEcal && dist > 0.1 ) {
          std::cout << "Warning - first ECAL cluster at a distance of " << dist << "from the closest track!" << std::endl;
          cout<<"\telement "<<elements[index]<<" linked with distance = "<< dist <<endl;
          cout<<"\telement "<<elements[iTrack]<<" linked with distance = "<< dist <<endl;
        }
        */ 

        // Check the presence of preshower clusters in the vicinity
        // Preshower cluster closer to another ECAL cluster are ignored.
        vector<double> ps1Ene(1,static_cast<double>(0.));
        associatePSClusters(index, reco::PFBlockElement::PS1, block, elements, linkData, active, ps1Ene);
        vector<double> ps2Ene(1,static_cast<double>(0.));
        associatePSClusters(index, reco::PFBlockElement::PS2, block, elements, linkData, active, ps2Ene);
        
        // Get the energy calibrated (for photons)
        bool crackCorrection = false;
        double ecalEnergy = calibration_->energyEm(*clusterRef,ps1Ene,ps2Ene,crackCorrection);
        if ( debug_ )
          std::cout << "Corrected ECAL(+PS) energy = " << ecalEnergy << std::endl;

        // Since the electrons were found beforehand, this track must be a hadron. Calibrate 
        // the energy under the hadron hypothesis.
        totalEcal += ecalEnergy;
        double previousCalibEcal = calibEcal;
        double previousSlopeEcal = slopeEcal;
        calibEcal = std::max(totalEcal,0.);
        calibHcal = 0.;
        calibration_->energyEmHad(trackMomentum,calibEcal,calibHcal,
                                  clusterRef->positionREP().Eta(),
                                  clusterRef->positionREP().Phi());
        if ( totalEcal > 0.) slopeEcal = calibEcal/totalEcal;

        if ( debug_ )
          std::cout << "The total calibrated energy so far amounts to = " << calibEcal << std::endl;
        

        // Stop the loop when adding more ECAL clusters ruins the compatibility
        if ( connectedToEcal && calibEcal - trackMomentum >= 0. ) {
        // if ( connectedToEcal && calibEcal - trackMomentum >=
        //     nSigmaECAL_*neutralHadronEnergyResolution(trackMomentum,clusterRef->positionREP().Eta())  ) { 
          calibEcal = previousCalibEcal;
          slopeEcal = previousSlopeEcal;
          totalEcal = calibEcal/slopeEcal;

          // Turn this last cluster in a photon 
          // (The PS clusters are already locked in "associatePSClusters")
          active[index] = false;

          // Find the associated tracks
          std::multimap<double, unsigned> assTracks;
          block.associatedElements( index,  linkData,
                                    assTracks,
                                    reco::PFBlockElement::TRACK,
                                    reco::PFBlock::LINKTEST_ALL );


          unsigned tmpe = reconstructCluster( *clusterRef, ecalEnergy ); 
          (*pfCandidates_)[tmpe].setEcalEnergy( clusterRef->energy(), ecalEnergy );
          (*pfCandidates_)[tmpe].setHcalEnergy( 0., 0. );
          (*pfCandidates_)[tmpe].setHoEnergy( 0., 0. );
          (*pfCandidates_)[tmpe].setPs1Energy( ps1Ene[0] );
          (*pfCandidates_)[tmpe].setPs2Energy( ps2Ene[0] );
          (*pfCandidates_)[tmpe].addElementInBlock( blockref, index );
          // Check that there is at least one track
          if(assTracks.size()) {
            (*pfCandidates_)[tmpe].addElementInBlock( blockref, assTracks.begin()->second );
            
            // Assign the position of the track at the ECAL entrance
            const math::XYZPointF& chargedPosition = 
              dynamic_cast<const reco::PFBlockElementTrack*>(&elements[assTracks.begin()->second])->positionAtECALEntrance();
            (*pfCandidates_)[tmpe].setPositionAtECALEntrance(chargedPosition);
          }
          break;
        }

        // Lock used clusters.
        connectedToEcal = true;
        iEcal = index;
        active[index] = false;
        for (unsigned ic=0; ic<tmpi.size();++ic)  
          (*pfCandidates_)[tmpi[ic]].addElementInBlock( blockref, iEcal ); 


      } // Loop ecal elements
    
      bool bNeutralProduced = false;
  
      // Create a photon if the ecal energy is too large
      if( connectedToEcal ) {
        reco::PFClusterRef pivotalRef = elements[iEcal].clusterRef();

        double neutralEnergy = calibEcal - trackMomentum;

        /*
        // Check if there are other tracks linked to that ECAL
        for(IE ie = sortedTracks.begin(); ie != sortedTracks.end() && neutralEnergy > 0; ++ie ) {
          unsigned jTrack = ie->second;
          
          // The loop is on the other tracks !
          if ( jTrack == iTrack ) continue;

          // Check if the ECAL closest to this track is the current ECAL
          // Otherwise ignore this track in the neutral energy determination
          std::multimap<double, unsigned> sortedECAL;
          block.associatedElements( jTrack,  linkData,
                                    sortedECAL,
                                    reco::PFBlockElement::ECAL,
                                    reco::PFBlock::LINKTEST_ALL );
          if ( sortedECAL.begin()->second != iEcal ) continue;
          
          // Check if this track is also linked to an HCAL 
          // (in which case it goes to the next loop and is ignored here)
          std::multimap<double, unsigned> sortedHCAL;
          block.associatedElements( jTrack,  linkData,
                                    sortedHCAL,
                                    reco::PFBlockElement::HCAL,
                                    reco::PFBlock::LINKTEST_ALL );
          if ( sortedHCAL.size() ) continue;
          
          // Otherwise, subtract this track momentum from the ECAL energy 
          reco::TrackRef trackRef = elements[jTrack].trackRef();
          neutralEnergy -= trackRef->p();

        } // End other tracks
        */
        
        // Add a photon is the energy excess is large enough
        double resol = neutralHadronEnergyResolution(trackMomentum,pivotalRef->positionREP().Eta());
        resol *= trackMomentum;
        if ( neutralEnergy > std::max(0.5,nSigmaECAL_*resol) ) {
          neutralEnergy /= slopeEcal;
          unsigned tmpj = reconstructCluster( *pivotalRef, neutralEnergy ); 
          (*pfCandidates_)[tmpj].setEcalEnergy( pivotalRef->energy(), neutralEnergy );
          (*pfCandidates_)[tmpj].setHcalEnergy( 0., 0. );
          (*pfCandidates_)[tmpj].setHoEnergy( 0., 0. );
          (*pfCandidates_)[tmpj].setPs1Energy( 0. );
          (*pfCandidates_)[tmpj].setPs2Energy( 0. );
          (*pfCandidates_)[tmpj].addElementInBlock(blockref, iEcal);
          bNeutralProduced = true;
          for (unsigned ic=0; ic<kTrack.size();++ic) 
            (*pfCandidates_)[tmpj].addElementInBlock( blockref, kTrack[ic] ); 
        } // End neutral energy

        // Set elements in blocks and ECAL energies to all tracks
        for (unsigned ic=0; ic<tmpi.size();++ic) { 
          
          // Skip muons
          if ( (*pfCandidates_)[tmpi[ic]].particleId() == reco::PFCandidate::mu ) continue; 

          double fraction = (*pfCandidates_)[tmpi[ic]].trackRef()->p()/trackMomentum;
          double ecalCal = bNeutralProduced ? 
            (calibEcal-neutralEnergy*slopeEcal)*fraction : calibEcal*fraction;
          double ecalRaw = totalEcal*fraction;

          if (debug_) cout << "The fraction after photon supression is " << fraction << " calibrated ecal = " << ecalCal << endl;

          (*pfCandidates_)[tmpi[ic]].setEcalEnergy( ecalRaw, ecalCal );
          (*pfCandidates_)[tmpi[ic]].setHcalEnergy( 0., 0. );
          (*pfCandidates_)[tmpi[ic]].setHoEnergy( 0., 0. );
          (*pfCandidates_)[tmpi[ic]].setPs1Energy( 0 );
          (*pfCandidates_)[tmpi[ic]].setPs2Energy( 0 );
          (*pfCandidates_)[tmpi[ic]].addElementInBlock( blockref, kTrack[ic] );
        }

      } // End connected ECAL

      // Fill the element_in_block for tracks that are eventually linked to no ECAL clusters at all.
      for (unsigned ic=0; ic<tmpi.size();++ic) { 
        const PFCandidate& pfc = (*pfCandidates_)[tmpi[ic]];
        const PFCandidate::ElementsInBlocks& eleInBlocks = pfc.elementsInBlocks();
        if ( eleInBlocks.size() == 0 ) { 
          if ( debug_ )std::cout << "Single track / Fill element in block! " << std::endl;
          (*pfCandidates_)[tmpi[ic]].addElementInBlock( blockref, kTrack[ic] );
        }
      }

    }


    // In case several HCAL elements are linked to this track, 
    // unlinking all of them except the closest. 
    for(IE ie = hcalElems.begin(); ie != hcalElems.end(); ++ie ) {
      
      unsigned index = ie->second;
      
      PFBlockElement::Type type = elements[index].type();

      if(debug_) {
        double dist = block.dist(iTrack,index,linkData,reco::PFBlock::LINKTEST_ALL);
        cout<<"\telement "<<elements[index]<<" linked with distance "<< dist <<endl;
      }

      assert( type == PFBlockElement::HCAL );
      
      // all hcal clusters except the closest 
      // will be unlinked from the track
      if( !hcalFound ) { // closest hcal
        if(debug_) 
          cout<<"\t\tclosest hcal cluster, doing nothing"<<endl;
        
        hcalFound = true;
        
        // active[index] = false;
        // hcalUsed.push_back( index );
      }
      else { // other associated hcal
        // unlink from the track
        if(debug_) 
          cout<<"\t\tsecondary hcal cluster. unlinking"<<endl;
        block.setLink( iTrack, index, -1., linkData,
                       PFBlock::LINKTEST_RECHIT );
      }
    } //loop hcal elements   
  } // end of loop 1 on elements iEle of any type



  // deal with HF. 
  if( !(hfEmIs.empty() &&  hfHadIs.empty() ) ) {
    // there is at least one HF element in this block. 
    // so all elements must be HF.
    assert( hfEmIs.size() + hfHadIs.size() == elements.size() );

    if( elements.size() == 1 ) {
      //Auguste: HAD-only calibration here
      reco::PFClusterRef clusterRef = elements[0].clusterRef();
      double energyHF = 0.;
      double uncalibratedenergyHF = 0.;
      unsigned tmpi = 0;
      switch( clusterRef->layer() ) {
      case PFLayer::HF_EM:
        // do EM-only calibration here
        energyHF = clusterRef->energy();
        uncalibratedenergyHF = energyHF;
        if(thepfEnergyCalibrationHF_->getcalibHF_use()==true){
          energyHF = thepfEnergyCalibrationHF_->energyEm(uncalibratedenergyHF,
                                                         clusterRef->positionREP().Eta(),
                                                         clusterRef->positionREP().Phi()); 
        }
        tmpi = reconstructCluster( *clusterRef, energyHF );     
        (*pfCandidates_)[tmpi].setEcalEnergy( uncalibratedenergyHF, energyHF );
        (*pfCandidates_)[tmpi].setHcalEnergy( 0., 0.);
        (*pfCandidates_)[tmpi].setHoEnergy( 0., 0.);
        (*pfCandidates_)[tmpi].setPs1Energy( 0. );
        (*pfCandidates_)[tmpi].setPs2Energy( 0. );
        (*pfCandidates_)[tmpi].addElementInBlock( blockref, hfEmIs[0] );
        //std::cout << "HF EM alone ! " << energyHF << std::endl;
        break;
      case PFLayer::HF_HAD:
        // do HAD-only calibration here
        energyHF = clusterRef->energy();
        uncalibratedenergyHF = energyHF;
        if(thepfEnergyCalibrationHF_->getcalibHF_use()==true){
          energyHF = thepfEnergyCalibrationHF_->energyHad(uncalibratedenergyHF,
                                                         clusterRef->positionREP().Eta(),
                                                         clusterRef->positionREP().Phi()); 
        }
        tmpi = reconstructCluster( *clusterRef, energyHF );     
        (*pfCandidates_)[tmpi].setHcalEnergy( uncalibratedenergyHF, energyHF );
        (*pfCandidates_)[tmpi].setEcalEnergy( 0., 0.);
        (*pfCandidates_)[tmpi].setHoEnergy( 0., 0.);
        (*pfCandidates_)[tmpi].setPs1Energy( 0. );
        (*pfCandidates_)[tmpi].setPs2Energy( 0. );
        (*pfCandidates_)[tmpi].addElementInBlock( blockref, hfHadIs[0] );
        //std::cout << "HF Had alone ! " << energyHF << std::endl;
        break;
      default:
        assert(0);
      }
    }
    else if(  elements.size() == 2 ) {
      //Auguste: EM + HAD calibration here
      reco::PFClusterRef c0 = elements[0].clusterRef();
      reco::PFClusterRef c1 = elements[1].clusterRef();
      // 2 HF elements. Must be in each layer.
      reco::PFClusterRef cem = c1;
      reco::PFClusterRef chad = c0;
      
      if( c0->layer()== PFLayer::HF_EM ) {
        if ( c1->layer()== PFLayer::HF_HAD ) {
          cem = c0; chad = c1;
        } else { 
          assert(0);
        }
        // do EM+HAD calibration here
        double energyHfEm = cem->energy();
        double energyHfHad = chad->energy();
        double uncalibratedenergyHFEm = energyHfEm;
        double uncalibratedenergyHFHad = energyHfHad;
        if(thepfEnergyCalibrationHF_->getcalibHF_use()==true){

          energyHfEm = thepfEnergyCalibrationHF_->energyEmHad(uncalibratedenergyHFEm,
                                                             0.0,
                                                             c0->positionREP().Eta(),
                                                             c0->positionREP().Phi()); 
          energyHfHad = thepfEnergyCalibrationHF_->energyEmHad(0.0,
                                                              uncalibratedenergyHFHad,
                                                              c1->positionREP().Eta(),
                                                              c1->positionREP().Phi()); 
        }
        unsigned tmpi = reconstructCluster( *chad, energyHfEm+energyHfHad );     
        (*pfCandidates_)[tmpi].setEcalEnergy( uncalibratedenergyHFEm, energyHfEm );
        (*pfCandidates_)[tmpi].setHcalEnergy( uncalibratedenergyHFHad, energyHfHad);
        (*pfCandidates_)[tmpi].setHoEnergy( 0., 0.);
        (*pfCandidates_)[tmpi].setPs1Energy( 0. );
        (*pfCandidates_)[tmpi].setPs2Energy( 0. );
        (*pfCandidates_)[tmpi].addElementInBlock( blockref, hfEmIs[0] );
        (*pfCandidates_)[tmpi].addElementInBlock( blockref, hfHadIs[0] );
        //std::cout << "HF EM+HAD found ! " << energyHfEm << " " << energyHfHad << std::endl;

      }
      else {
        cerr<<"Warning: 2 elements, but not 1 HFEM and 1 HFHAD"<<endl;
        cerr<<block<<endl;
//      assert( c1->layer()== PFLayer::HF_EM &&
//              c0->layer()== PFLayer::HF_HAD );
      }      
    }
    else {
      // 1 HF element in the block, 
      // but number of elements not equal to 1 or 2 
      cerr<<"Warning: HF, but n elem different from 1 or 2"<<endl;
      cerr<<block<<endl;
//       assert(0);
//       cerr<<"not ready for navigation in the HF!"<<endl;
    }
  }



  if(debug_) { 
    cout<<endl;
    cout<<endl<<"--------------- loop hcal ---------------------"<<endl;
  }
  
  // track information:
  // trackRef
  // ecal energy
  // hcal energy
  // rescale 

  for(unsigned i=0; i<hcalIs.size(); i++) {
    
    unsigned iHcal= hcalIs[i];    
    PFBlockElement::Type type = elements[iHcal].type();

    assert( type == PFBlockElement::HCAL );

    if(debug_) cout<<endl<<elements[iHcal]<<endl;

    // vector<unsigned> elementIndices;
    // elementIndices.push_back(iHcal);

    // associated tracks
    std::multimap<double, unsigned> sortedTracks;
    block.associatedElements( iHcal,  linkData,
                              sortedTracks,
                              reco::PFBlockElement::TRACK,
                              reco::PFBlock::LINKTEST_ALL );

    std::multimap< unsigned, std::pair<double, unsigned> > associatedEcals;

    std::map< unsigned, std::pair<double, double> > associatedPSs;
 
    std::multimap<double, std::pair<unsigned,bool> > associatedTracks;
    
    // A temporary maps for ECAL satellite clusters
    std::multimap<double,std::pair<unsigned,math::XYZVector> > ecalSatellites;
    std::pair<unsigned,math::XYZVector> fakeSatellite = make_pair(iHcal,math::XYZVector(0.,0.,0.));
    ecalSatellites.insert( make_pair(-1., fakeSatellite) );

    std::multimap< unsigned, std::pair<double, unsigned> > associatedHOs;
    
    PFClusterRef hclusterref = elements[iHcal].clusterRef();
    assert(!hclusterref.isNull() );


    //if there is no track attached to that HCAL, then do not
    //reconstruct an HCAL alone, check if it can be recovered 
    //first
    
    // if no associated tracks, create a neutral hadron
    //if(sortedTracks.empty() ) {

    if( sortedTracks.empty() ) {
      if(debug_) 
        cout<<"\tno associated tracks, keep for later"<<endl;
      continue;
    }

    // Lock the HCAL cluster
    active[iHcal] = false;


    // in the following, tracks are associated to this hcal cluster.
    // will look for an excess of energy in the calorimeters w/r to 
    // the charged energy, and turn this excess into a neutral hadron or 
    // a photon.

    if(debug_) cout<<"\t"<<sortedTracks.size()<<" associated tracks:"<<endl;

    double totalChargedMomentum = 0;
    double sumpError2 = 0.;
    double totalHO = 0.;
    double totalEcal = 0.;
    double totalHcal = hclusterref->energy();
    vector<double> hcalP;
    vector<double> hcalDP;
    vector<unsigned> tkIs;
    double maxDPovP = -9999.;
    
    //Keep track of how much energy is assigned to calorimeter-vs-track energy/momentum excess
    vector< unsigned > chargedHadronsIndices;
    vector< unsigned > chargedHadronsInBlock;
    double mergedNeutralHadronEnergy = 0;
    double mergedPhotonEnergy = 0;
    double muonHCALEnergy = 0.;
    double muonECALEnergy = 0.;
    double muonHCALError = 0.;
    double muonECALError = 0.;
    unsigned nMuons = 0; 


    math::XYZVector photonAtECAL(0.,0.,0.);
    math::XYZVector hadronDirection(hclusterref->position().X(),
                                    hclusterref->position().Y(),
                                    hclusterref->position().Z());
    hadronDirection = hadronDirection.Unit();
    math::XYZVector hadronAtECAL = totalHcal * hadronDirection;

    // Loop over all tracks associated to this HCAL cluster
    for(IE ie = sortedTracks.begin(); ie != sortedTracks.end(); ++ie ) {

      unsigned iTrack = ie->second;

      // Check the track has not already been used (e.g., in electrons, conversions...)
      if ( !active[iTrack] ) continue;
      // Sanity check 1
      PFBlockElement::Type type = elements[iTrack].type();
      assert( type == reco::PFBlockElement::TRACK );
      // Sanity check 2
      reco::TrackRef trackRef = elements[iTrack].trackRef();
      assert( !trackRef.isNull() ); 

      // The direction at ECAL entrance
      const math::XYZPointF& chargedPosition = 
        dynamic_cast<const reco::PFBlockElementTrack*>(&elements[iTrack])->positionAtECALEntrance();
      math::XYZVector chargedDirection(chargedPosition.X(),chargedPosition.Y(),chargedPosition.Z());
      chargedDirection = chargedDirection.Unit();

      // look for ECAL elements associated to iTrack (associated to iHcal)
      std::multimap<double, unsigned> sortedEcals;
      block.associatedElements( iTrack,  linkData,
                                sortedEcals,
                                reco::PFBlockElement::ECAL,
                                reco::PFBlock::LINKTEST_ALL );

      if(debug_) cout<<"\t\t\tnumber of Ecal elements linked to this track: "
                     <<sortedEcals.size()<<endl;     

      // look for HO elements associated to iTrack (associated to iHcal)
      std::multimap<double, unsigned> sortedHOs;
      if (useHO_) {
        block.associatedElements( iTrack,  linkData,
                                  sortedHOs,
                                  reco::PFBlockElement::HO,
                                  reco::PFBlock::LINKTEST_ALL );

      }
      if(debug_) 
        cout<<"PFAlgo : number of HO elements linked to this track: "
            <<sortedHOs.size()<<endl;
      
      // Create a PF Candidate right away if the track is a tight muon
      reco::MuonRef muonRef = elements[iTrack].muonRef();
      bool thisIsAMuon = PFMuonAlgo::isMuon(elements[iTrack]);
      bool thisIsAnIsolatedMuon = PFMuonAlgo::isIsolatedMuon(elements[iTrack]);
      bool thisIsALooseMuon = false;

      if(!thisIsAMuon ) {
        thisIsALooseMuon = PFMuonAlgo::isLooseMuon(elements[iTrack]);
      }
      if ( thisIsAMuon ) {
        if ( debug_ ) { 
          std::cout << "\t\tThis track is identified as a muon - remove it from the stack" << std::endl;
          std::cout << "\t\t" << elements[iTrack] << std::endl;
        }

        // Create a muon.
        unsigned tmpi = reconstructTrack( elements[iTrack] );
        (*pfCandidates_)[tmpi].addElementInBlock( blockref, iTrack );
        (*pfCandidates_)[tmpi].addElementInBlock( blockref, iHcal );
        double muonHcal = std::min(muonHCAL_[0]+muonHCAL_[1],totalHcal);

        // if muon is isolated and muon momentum exceeds the calo energy, absorb the calo energy        
        // rationale : there has been a EM showering by the muon in the calorimeter - or the coil -
        // and we don't want to double count the energy 
        bool letMuonEatCaloEnergy = false;

        if(thisIsAnIsolatedMuon){
          // The factor 1.3 is the e/pi factor in HCAL...
          double totalCaloEnergy = totalHcal / 1.30;
          unsigned iEcal = 0;
          if( !sortedEcals.empty() ) { 
            iEcal = sortedEcals.begin()->second; 
            PFClusterRef eclusterref = elements[iEcal].clusterRef();
            totalCaloEnergy += eclusterref->energy();
          }

          if (useHO_) {
            // The factor 1.3 is assumed to be the e/pi factor in HO, too.
            unsigned iHO = 0;
            if( !sortedHOs.empty() ) { 
              iHO = sortedHOs.begin()->second; 
              PFClusterRef eclusterref = elements[iHO].clusterRef();
              totalCaloEnergy += eclusterref->energy() / 1.30;
            }
          }

          // std::cout << "muon p / total calo = " << muonRef->p() << " "  << (pfCandidates_->back()).p() << " " << totalCaloEnergy << std::endl;
          //if(muonRef->p() > totalCaloEnergy ) letMuonEatCaloEnergy = true;
          if( (pfCandidates_->back()).p() > totalCaloEnergy ) letMuonEatCaloEnergy = true;
        }

        if(letMuonEatCaloEnergy) muonHcal = totalHcal;
        double muonEcal =0.;
        unsigned iEcal = 0;
        if( !sortedEcals.empty() ) { 
          iEcal = sortedEcals.begin()->second; 
          PFClusterRef eclusterref = elements[iEcal].clusterRef();
          (*pfCandidates_)[tmpi].addElementInBlock( blockref, iEcal);
          muonEcal = std::min(muonECAL_[0]+muonECAL_[1],eclusterref->energy());
          if(letMuonEatCaloEnergy) muonEcal = eclusterref->energy();
          // If the muon expected energy accounts for the whole ecal cluster energy, lock the ecal cluster
          if ( eclusterref->energy() - muonEcal  < 0.2 ) active[iEcal] = false;
          (*pfCandidates_)[tmpi].setEcalEnergy(eclusterref->energy(), muonEcal);
        }
        unsigned iHO = 0;
        double muonHO =0.;
        if (useHO_) {
          if( !sortedHOs.empty() ) { 
            iHO = sortedHOs.begin()->second; 
            PFClusterRef hoclusterref = elements[iHO].clusterRef();
            (*pfCandidates_)[tmpi].addElementInBlock( blockref, iHO);
            muonHO = std::min(muonHO_[0]+muonHO_[1],hoclusterref->energy());
            if(letMuonEatCaloEnergy) muonHO = hoclusterref->energy();
            // If the muon expected energy accounts for the whole HO cluster energy, lock the HO cluster
            if ( hoclusterref->energy() - muonHO  < 0.2 ) active[iHO] = false;      
            (*pfCandidates_)[tmpi].setHcalEnergy(totalHcal, muonHcal);
            (*pfCandidates_)[tmpi].setHoEnergy(hoclusterref->energy(), muonHO);
          }
        } else {
          (*pfCandidates_)[tmpi].setHcalEnergy(totalHcal, muonHcal);
        }

        if(letMuonEatCaloEnergy){
          muonHCALEnergy += totalHcal;
          if (useHO_) muonHCALEnergy +=muonHO;
          muonHCALError += 0.;
          muonECALEnergy += muonEcal;
          muonECALError += 0.;
          photonAtECAL -= muonEcal*chargedDirection;
          hadronAtECAL -= totalHcal*chargedDirection;
          if ( !sortedEcals.empty() ) active[iEcal] = false;
          active[iHcal] = false;
          if (useHO_ && !sortedHOs.empty() ) active[iHO] = false;
        }
        else{
        // Estimate of the energy deposit & resolution in the calorimeters
          muonHCALEnergy += muonHCAL_[0];
          muonHCALError += muonHCAL_[1]*muonHCAL_[1];
          if ( muonHO > 0. ) { 
            muonHCALEnergy += muonHO_[0];
            muonHCALError += muonHO_[1]*muonHO_[1];
          }
          muonECALEnergy += muonECAL_[0];
          muonECALError += muonECAL_[1]*muonECAL_[1];
          // ... as well as the equivalent "momentum" at ECAL entrance
          photonAtECAL -= muonECAL_[0]*chargedDirection;
          hadronAtECAL -= muonHCAL_[0]*chargedDirection;
        }

        // Remove it from the stack
        active[iTrack] = false;
        // Go to next track
        continue;
      }
 
      //

      if(debug_) cout<<"\t\t"<<elements[iTrack]<<endl;

      // introduce  tracking errors 
      double trackMomentum = trackRef->p();
      totalChargedMomentum += trackMomentum; 

      // If the track is not a tight muon, but still resembles a muon
      // keep it for later for blocks with too large a charged energy
      if ( thisIsALooseMuon && !thisIsAMuon ) nMuons += 1;      

      // ... and keep anyway the pt error for possible fake rejection
      // ... blow up errors of 5th anf 4th iteration, to reject those
      // ... tracks first (in case it's needed)
      double Dpt = trackRef->ptError();
      double blowError = 1.;
      switch (trackRef->algo()) {
      case TrackBase::ctf:
      case TrackBase::iter0:
      case TrackBase::iter1:
      case TrackBase::iter2:
      case TrackBase::iter3:
      case TrackBase::iter4:
        blowError = 1.;
        break;
      case TrackBase::iter5:
        blowError = factors45_[0];
        break;
      case TrackBase::iter6:
        blowError = factors45_[1];
        break;
      default:
        blowError = 1E9;
        break;
      }
      // except if it is from an interaction
      bool isPrimaryOrSecondary = isFromSecInt(elements[iTrack], "all");

      if ( isPrimaryOrSecondary ) blowError = 1.;

      std::pair<unsigned,bool> tkmuon = make_pair(iTrack,thisIsALooseMuon);
      associatedTracks.insert( make_pair(-Dpt*blowError, tkmuon) );     
 
      // Also keep the total track momentum error for comparison with the calo energy
      double Dp = trackRef->qoverpError()*trackMomentum*trackMomentum;
      sumpError2 += Dp*Dp;

      bool connectedToEcal = false; // Will become true if there is at least one ECAL cluster connected
      if( !sortedEcals.empty() ) 
        { // start case: at least one ecal element associated to iTrack
          
          // Loop over all connected ECAL clusters
          for ( IE iec=sortedEcals.begin(); 
                iec!=sortedEcals.end(); ++iec ) {

            unsigned iEcal = iec->second;
            double dist = iec->first;
            
            // Ignore ECAL cluters already used
            if( !active[iEcal] ) { 
              if(debug_) cout<<"\t\tcluster locked"<<endl;          
              continue;
            }

            // Sanity checks
            PFBlockElement::Type type = elements[ iEcal ].type();
            assert( type == PFBlockElement::ECAL );
            PFClusterRef eclusterref = elements[iEcal].clusterRef();
            assert(!eclusterref.isNull() );
              
            // Check if this ECAL is not closer to another track - ignore it in that case
            std::multimap<double, unsigned> sortedTracksEcal;
            block.associatedElements( iEcal,  linkData,
                                      sortedTracksEcal,
                                      reco::PFBlockElement::TRACK,
                                      reco::PFBlock::LINKTEST_ALL );
            unsigned jTrack = sortedTracksEcal.begin()->second;
            if ( jTrack != iTrack ) continue; 

            // double chi2Ecal = block.chi2(jTrack,iEcal,linkData,
            //                              reco::PFBlock::LINKTEST_ALL);
            double distEcal = block.dist(jTrack,iEcal,linkData,
                                         reco::PFBlock::LINKTEST_ALL);
                                         
            // totalEcal += eclusterref->energy();
            // float ecalEnergyUnCalibrated = eclusterref->energy();
            //std::cout << "Ecal Uncalibrated " << ecalEnergyUnCalibrated << std::endl;
            
            // check the presence of preshower clusters in the vicinity
            vector<double> ps1Ene(1,static_cast<double>(0.));
            associatePSClusters(iEcal, reco::PFBlockElement::PS1, 
                                block, elements, linkData, active, 
                                ps1Ene);
            vector<double> ps2Ene(1,static_cast<double>(0.));
            associatePSClusters(iEcal, reco::PFBlockElement::PS2, 
                                block, elements, linkData, active, 
                                ps2Ene); 
            std::pair<double,double> psEnes = make_pair(ps1Ene[0],
                                                        ps2Ene[0]);
            associatedPSs.insert(make_pair(iEcal,psEnes));
              
            // Calibrate the ECAL energy for photons
            bool crackCorrection = false;
            float ecalEnergyCalibrated = calibration_->energyEm(*eclusterref,ps1Ene,ps2Ene,crackCorrection);
            math::XYZVector photonDirection(eclusterref->position().X(),
                                            eclusterref->position().Y(),
                                            eclusterref->position().Z());
            photonDirection = photonDirection.Unit();

            if ( !connectedToEcal ) { // This is the closest ECAL cluster - will add its energy later
              
              if(debug_) cout<<"\t\t\tclosest: "
                             <<elements[iEcal]<<endl;
              
              connectedToEcal = true;
              // PJ 1st-April-09 : To be done somewhere !!! (Had to comment it, but it is needed)
              // currentChargedHadron.addElementInBlock( blockref, iEcal );
       
              std::pair<unsigned,math::XYZVector> satellite = 
                make_pair(iEcal,ecalEnergyCalibrated*photonDirection);
              ecalSatellites.insert( make_pair(-1., satellite) );

            } else { // Keep satellite clusters for later
              
              std::pair<unsigned,math::XYZVector> satellite = 
                make_pair(iEcal,ecalEnergyCalibrated*photonDirection);
              ecalSatellites.insert( make_pair(dist, satellite) );
              
            }
            
            std::pair<double, unsigned> associatedEcal 
              = make_pair( distEcal, iEcal );
            associatedEcals.insert( make_pair(iTrack, associatedEcal) );

          } // End loop ecal associated to iTrack
        } // end case: at least one ecal element associated to iTrack

      if( useHO_ && !sortedHOs.empty() ) 
        { // start case: at least one ho element associated to iTrack
          
          // Loop over all connected HO clusters
          for ( IE ieho=sortedHOs.begin(); ieho!=sortedHOs.end(); ++ieho ) {

            unsigned iHO = ieho->second;
            double distHO = ieho->first;
            
            // Ignore HO cluters already used
            if( !active[iHO] ) { 
              if(debug_) cout<<"\t\tcluster locked"<<endl;          
              continue;
            }
            
            // Sanity checks
            PFBlockElement::Type type = elements[ iHO ].type();
            assert( type == PFBlockElement::HO );
            PFClusterRef hoclusterref = elements[iHO].clusterRef();
            assert(!hoclusterref.isNull() );
            
            // Check if this HO is not closer to another track - ignore it in that case
            std::multimap<double, unsigned> sortedTracksHO;
            block.associatedElements( iHO,  linkData,
                                      sortedTracksHO,
                                      reco::PFBlockElement::TRACK,
                                      reco::PFBlock::LINKTEST_ALL );
            unsigned jTrack = sortedTracksHO.begin()->second;
            if ( jTrack != iTrack ) continue; 

            // double chi2HO = block.chi2(jTrack,iHO,linkData,
            //                              reco::PFBlock::LINKTEST_ALL);
            //double distHO = block.dist(jTrack,iHO,linkData,
            //                 reco::PFBlock::LINKTEST_ALL);
            
            // Increment the total energy by the energy of this HO cluster
            totalHO += hoclusterref->energy();
            active[iHO] = false;
            // Keep track for later reference in the PFCandidate.
            std::pair<double, unsigned> associatedHO
              = make_pair( distHO, iHO );
            associatedHOs.insert( make_pair(iTrack, associatedHO) );
            
          } // End loop ho associated to iTrack
        } // end case: at least one ho element associated to iTrack
      
      
    } // end loop on tracks associated to hcal element iHcal

    // Include totalHO in totalHCAL for the time being (it will be calibrated as HCAL energy)
    totalHcal += totalHO;
    
    // test compatibility between calo and tracker. //////////////

    double caloEnergy = 0.;
    double slopeEcal = 1.0;
    double calibEcal = 0.;
    double calibHcal = 0.;
    hadronDirection = hadronAtECAL.Unit();

    // Determine the expected calo resolution from the total charged momentum
    double Caloresolution = neutralHadronEnergyResolution( totalChargedMomentum, hclusterref->positionREP().Eta());    
    Caloresolution *= totalChargedMomentum;
    // Account for muons
    Caloresolution = std::sqrt(Caloresolution*Caloresolution + muonHCALError + muonECALError);
    totalEcal -= std::min(totalEcal,muonECALEnergy);
    totalHcal -= std::min(totalHcal,muonHCALEnergy);
    if ( totalEcal < 1E-9 ) photonAtECAL = math::XYZVector(0.,0.,0.);
    if ( totalHcal < 1E-9 ) hadronAtECAL = math::XYZVector(0.,0.,0.);

    // Loop over all ECAL satellites, starting for the closest to the various tracks
    // and adding other satellites until saturation of the total track momentum
    // Note : for code simplicity, the first element of the loop is the HCAL cluster
    // with 0 energy in the ECAL
    for ( IS is = ecalSatellites.begin(); is != ecalSatellites.end(); ++is ) { 

      // Add the energy of this ECAL cluster
      double previousCalibEcal = calibEcal;
      double previousCalibHcal = calibHcal;
      double previousCaloEnergy = caloEnergy;
      double previousSlopeEcal = slopeEcal;
      math::XYZVector previousHadronAtECAL = hadronAtECAL;
      //
      totalEcal += sqrt(is->second.second.Mag2());
      photonAtECAL += is->second.second;
      calibEcal = std::max(0.,totalEcal);
      calibHcal = std::max(0.,totalHcal);
      hadronAtECAL = calibHcal * hadronDirection;
      // Calibrate ECAL and HCAL energy under the hadron hypothesis.
      calibration_->energyEmHad(totalChargedMomentum,calibEcal,calibHcal,
                                hclusterref->positionREP().Eta(),
                                hclusterref->positionREP().Phi());
      caloEnergy = calibEcal+calibHcal;
      if ( totalEcal > 0.) slopeEcal = calibEcal/totalEcal;

      hadronAtECAL = calibHcal * hadronDirection; 
      
      // Continue looping until all closest clusters are exhausted and as long as
      // the calorimetric energy does not saturate the total momentum.
      if ( is->first < 0. || caloEnergy - totalChargedMomentum <= 0. ) { 
        if(debug_) cout<<"\t\t\tactive, adding "<<is->second.second
                       <<" to ECAL energy, and locking"<<endl;
        active[is->second.first] = false;
        continue;
      }

      // Otherwise, do not consider the last cluster examined and exit.
      // active[is->second.first] = true;
      totalEcal -= sqrt(is->second.second.Mag2());
      photonAtECAL -= is->second.second;
      calibEcal = previousCalibEcal;
      calibHcal = previousCalibHcal;
      hadronAtECAL = previousHadronAtECAL;
      slopeEcal = previousSlopeEcal;
      caloEnergy = previousCaloEnergy;

      break;

    }

    // Sanity check !
    assert(caloEnergy>=0);


    // And now check for hadronic energy excess...

    //colin: resolution should be measured on the ecal+hcal case. 
    // however, the result will be close. 
    // double Caloresolution = neutralHadronEnergyResolution( caloEnergy );
    // Caloresolution *= caloEnergy;
    // PJ The resolution is on the expected charged calo energy !
    //double Caloresolution = neutralHadronEnergyResolution( totalChargedMomentum, hclusterref->positionREP().Eta());
    //Caloresolution *= totalChargedMomentum;
    // that of the charged particles linked to the cluster!
    double TotalError = sqrt(sumpError2 + Caloresolution*Caloresolution);

    /* */
    if ( totalChargedMomentum - caloEnergy > nSigmaTRACK_*Caloresolution ) {    

      /* 
      cout<<"\tCompare Calo Energy to total charged momentum "<<endl;
      cout<<"\t\tsum p    = "<<totalChargedMomentum<<" +- "<<sqrt(sumpError2)<<endl;
      cout<<"\t\tsum ecal = "<<totalEcal<<endl;
      cout<<"\t\tsum hcal = "<<totalHcal<<endl;
      cout<<"\t\t => Calo Energy = "<<caloEnergy<<" +- "<<Caloresolution<<endl;
      cout<<"\t\t => Calo Energy- total charged momentum = "
         <<caloEnergy-totalChargedMomentum<<" +- "<<TotalError<<endl;
      cout<<endl;
      cout << "\t\tNumber/momentum of muons found in the block : " << nMuons << std::endl;
      */
      // First consider loose muons
      if ( nMuons > 0 ) { 
        for ( IT it = associatedTracks.begin(); it != associatedTracks.end(); ++it ) { 
          unsigned iTrack = it->second.first;
          // Only active tracks
          if ( !active[iTrack] ) continue;
          // Only muons
          if ( !it->second.second ) continue;

          bool isTrack = elements[it->second.first].muonRef()->isTrackerMuon();
          double trackMomentum = elements[it->second.first].trackRef()->p();

          // look for ECAL elements associated to iTrack (associated to iHcal)
          std::multimap<double, unsigned> sortedEcals;
          block.associatedElements( iTrack,  linkData,
                                    sortedEcals,
                                    reco::PFBlockElement::ECAL,
                                    reco::PFBlock::LINKTEST_ALL );
          std::multimap<double, unsigned> sortedHOs;
          block.associatedElements( iTrack,  linkData,
                                    sortedHOs,
                                    reco::PFBlockElement::HO,
                                    reco::PFBlock::LINKTEST_ALL );
          // Add the muon to the PF candidate list
          unsigned tmpi = reconstructTrack( elements[iTrack] );
          (*pfCandidates_)[tmpi].addElementInBlock( blockref, iTrack );
          (*pfCandidates_)[tmpi].addElementInBlock( blockref, iHcal );
          double muonHcal = std::min(muonHCAL_[0]+muonHCAL_[1],totalHcal-totalHO);
          double muonHO = 0.;
          (*pfCandidates_)[tmpi].setHcalEnergy(totalHcal,muonHcal);
          if( !sortedEcals.empty() ) { 
            unsigned iEcal = sortedEcals.begin()->second; 
            PFClusterRef eclusterref = elements[iEcal].clusterRef();
            (*pfCandidates_)[tmpi].addElementInBlock( blockref, iEcal);
            double muonEcal = std::min(muonECAL_[0]+muonECAL_[1],eclusterref->energy());
            (*pfCandidates_)[tmpi].setEcalEnergy(eclusterref->energy(),muonEcal);
          }
          if( useHO_ && !sortedHOs.empty() ) { 
            unsigned iHO = sortedHOs.begin()->second; 
            PFClusterRef hoclusterref = elements[iHO].clusterRef();
            (*pfCandidates_)[tmpi].addElementInBlock( blockref, iHO);
            muonHO = std::min(muonHO_[0]+muonHO_[1],hoclusterref->energy());
            (*pfCandidates_)[tmpi].setHcalEnergy(totalHcal-totalHO,muonHcal);
            (*pfCandidates_)[tmpi].setHoEnergy(hoclusterref->energy(),muonHO);
          }
          reco::PFCandidate::ParticleType particleType = reco::PFCandidate::mu;
          (*pfCandidates_)[tmpi].setParticleType(particleType);

          if(PFMuonAlgo::isGlobalLooseMuon(elements[it->second.first].muonRef())){
            // Take the best pt measurement
            reco::TrackRef bestMuTrack =  elements[it->second.first].muonRef()->muonBestTrack();
            if(!useBestMuonTrack_)
              bestMuTrack = elements[it->second.first].muonRef()->combinedMuon();

            double muonMomentum = bestMuTrack->p();
            double muonPtError  = bestMuTrack->ptError();

            double staMomentum = elements[it->second.first].muonRef()->standAloneMuon()->p();
            double staPtError = elements[it->second.first].muonRef()->standAloneMuon()->ptError();
            double trackPtError = elements[it->second.first].trackRef()->ptError();


            double globalCorr = 1;

            // for loose muons we choose whichever fit has the best error, since the track is often mis-recontructed or there are problems with arbitration
            if(!isTrack || muonPtError < trackPtError || staPtError < trackPtError){          
              if(muonPtError < staPtError) globalCorr = muonMomentum/trackMomentum;
              else globalCorr = staMomentum/trackMomentum;
            }
            (*pfCandidates_)[tmpi].rescaleMomentum(globalCorr);
            
            if (debug_){
              std::cout << "\tElement  " << elements[iTrack] << std::endl 
                        << "PFAlgo: particle type set to muon (global, loose), pT = "<<elements[it->second.first].muonRef()->pt()<<std::endl; 
              PFMuonAlgo::printMuonProperties(elements[it->second.first].muonRef());
            }
          }
          else{
            if (debug_){
              std::cout << "\tElement  " << elements[iTrack] << std::endl 
                        << "PFAlgo: particle type set to muon (tracker, loose), pT = "<<elements[it->second.first].muonRef()->pt()<<std::endl;
              PFMuonAlgo::printMuonProperties(elements[it->second.first].muonRef()); 
            }
          }
          
          // Remove it from the block
          const math::XYZPointF& chargedPosition = 
            dynamic_cast<const reco::PFBlockElementTrack*>(&elements[it->second.first])->positionAtECALEntrance();        
          math::XYZVector chargedDirection(chargedPosition.X(), chargedPosition.Y(), chargedPosition.Z());
          chargedDirection = chargedDirection.Unit();
          totalChargedMomentum -= trackMomentum;
          // Update the calo energies
          if ( totalEcal > 0. ) 
            calibEcal -= std::min(calibEcal,muonECAL_[0]*calibEcal/totalEcal);
          if ( totalHcal > 0. ) 
            calibHcal -= std::min(calibHcal,muonHCAL_[0]*calibHcal/totalHcal);
          totalEcal -= std::min(totalEcal,muonECAL_[0]);
          totalHcal -= std::min(totalHcal,muonHCAL_[0]);
          if ( totalEcal > muonECAL_[0] ) photonAtECAL -= muonECAL_[0] * chargedDirection;
          if ( totalHcal > muonHCAL_[0] ) hadronAtECAL -= muonHCAL_[0]*calibHcal/totalHcal * chargedDirection;
          caloEnergy = calibEcal+calibHcal;
          muonHCALEnergy += muonHCAL_[0];
          muonHCALError += muonHCAL_[1]*muonHCAL_[1];
          if ( muonHO > 0. ) { 
            muonHCALEnergy += muonHO_[0];
            muonHCALError += muonHO_[1]*muonHO_[1];
            if ( totalHcal > 0. ) { 
              calibHcal -= std::min(calibHcal,muonHO_[0]*calibHcal/totalHcal);
              totalHcal -= std::min(totalHcal,muonHO_[0]);
            }
          }
          muonECALEnergy += muonECAL_[0];
          muonECALError += muonECAL_[1]*muonECAL_[1];
          active[iTrack] = false;
          // Stop the loop whenever enough muons are removed
          //Commented out: Keep looking for muons since they often come in pairs -Matt
          //if ( totalChargedMomentum < caloEnergy ) break;     
        }
        // New calo resolution.
        Caloresolution = neutralHadronEnergyResolution( totalChargedMomentum, hclusterref->positionREP().Eta());    
        Caloresolution *= totalChargedMomentum;
        Caloresolution = std::sqrt(Caloresolution*Caloresolution + muonHCALError + muonECALError);
      }
    }
    /*    
    if(debug_){
      cout<<"\tBefore Cleaning "<<endl;
      cout<<"\tCompare Calo Energy to total charged momentum "<<endl;
      cout<<"\t\tsum p    = "<<totalChargedMomentum<<" +- "<<sqrt(sumpError2)<<endl;
      cout<<"\t\tsum ecal = "<<totalEcal<<endl;
      cout<<"\t\tsum hcal = "<<totalHcal<<endl;
      cout<<"\t\t => Calo Energy = "<<caloEnergy<<" +- "<<Caloresolution<<endl;
      cout<<"\t\t => Calo Energy- total charged momentum = "
          <<caloEnergy-totalChargedMomentum<<" +- "<<TotalError<<endl;
      cout<<endl;
    }
    */

    // Second consider bad tracks (if still needed after muon removal)
    unsigned corrTrack = 10000000;
    double corrFact = 1.;

    if (rejectTracks_Bad_ && 
        totalChargedMomentum - caloEnergy > nSigmaTRACK_*Caloresolution) { 
      
      for ( IT it = associatedTracks.begin(); it != associatedTracks.end(); ++it ) { 
        unsigned iTrack = it->second.first;
        // Only active tracks
        if ( !active[iTrack] ) continue;
        const reco::TrackRef& trackref = elements[it->second.first].trackRef();

        double dptRel = fabs(it->first)/trackref->pt()*100;
        bool isSecondary = isFromSecInt(elements[iTrack], "secondary");
        bool isPrimary = isFromSecInt(elements[iTrack], "primary");

        if ( isSecondary && dptRel < dptRel_DispVtx_ ) continue;
        // Consider only bad tracks
        if ( fabs(it->first) < ptError_ ) break;
        // What would become the block charged momentum if this track were removed
        double wouldBeTotalChargedMomentum = 
          totalChargedMomentum - trackref->p();
        // Reject worst tracks, as long as the total charged momentum 
        // is larger than the calo energy

        if ( wouldBeTotalChargedMomentum > caloEnergy ) { 

          if (debug_ && isSecondary) {
            cout << "In bad track rejection step dptRel = " << dptRel << " dptRel_DispVtx_ = " << dptRel_DispVtx_ << endl;
            cout << "The calo energy would be still smaller even without this track but it is attached to a NI"<< endl;
          }


          if(isPrimary || (isSecondary && dptRel < dptRel_DispVtx_)) continue;
          active[iTrack] = false;
          totalChargedMomentum = wouldBeTotalChargedMomentum;
          if ( debug_ )
            std::cout << "\tElement  " << elements[iTrack] 
                      << " rejected (Dpt = " << -it->first 
                      << " GeV/c, algo = " << trackref->algo() << ")" << std::endl;
        // Just rescale the nth worst track momentum to equalize the calo energy
        } else {        
          if(isPrimary) break;
          corrTrack = iTrack;
          corrFact = (caloEnergy - wouldBeTotalChargedMomentum)/elements[it->second.first].trackRef()->p();
          if ( trackref->p()*corrFact < 0.05 ) { 
            corrFact = 0.;
            active[iTrack] = false;
          }
          totalChargedMomentum -= trackref->p()*(1.-corrFact);
          if ( debug_ ) 
            std::cout << "\tElement  " << elements[iTrack] 
                      << " (Dpt = " << -it->first 
                      << " GeV/c, algo = " << trackref->algo() 
                      << ") rescaled by " << corrFact 
                      << " Now the total charged momentum is " << totalChargedMomentum  << endl;
          break;
        }
      }
    }

    // New determination of the calo and track resolution avec track deletion/rescaling.
    Caloresolution = neutralHadronEnergyResolution( totalChargedMomentum, hclusterref->positionREP().Eta());    
    Caloresolution *= totalChargedMomentum;
    Caloresolution = std::sqrt(Caloresolution*Caloresolution + muonHCALError + muonECALError);

    // Check if the charged momentum is still very inconsistent with the calo measurement.
    // In this case, just drop all tracks from 4th and 5th iteration linked to this block
    
    if ( rejectTracks_Step45_ &&
         sortedTracks.size() > 1 && 
         totalChargedMomentum - caloEnergy > nSigmaTRACK_*Caloresolution ) { 
      for ( IT it = associatedTracks.begin(); it != associatedTracks.end(); ++it ) { 
        unsigned iTrack = it->second.first;
        reco::TrackRef trackref = elements[iTrack].trackRef();
        if ( !active[iTrack] ) continue;

        double dptRel = fabs(it->first)/trackref->pt()*100;
        bool isPrimaryOrSecondary = isFromSecInt(elements[iTrack], "all");




        if (isPrimaryOrSecondary && dptRel < dptRel_DispVtx_) continue;
        //
        switch (trackref->algo()) {
        case TrackBase::ctf:
        case TrackBase::iter0:
        case TrackBase::iter1:
        case TrackBase::iter2:
        case TrackBase::iter3:
        case TrackBase::iter4:
          break;
        case TrackBase::iter5:
        case TrackBase::iter6:
          active[iTrack] = false;       
          totalChargedMomentum -= trackref->p();
          
          if ( debug_ ) 
            std::cout << "\tElement  " << elements[iTrack] 
                      << " rejected (Dpt = " << -it->first 
                      << " GeV/c, algo = " << trackref->algo() << ")" << std::endl;
          break;
        default:
          break;
        }
      }
    }

    // New determination of the calo and track resolution avec track deletion/rescaling.
    Caloresolution = neutralHadronEnergyResolution( totalChargedMomentum, hclusterref->positionREP().Eta());    
    Caloresolution *= totalChargedMomentum;
    Caloresolution = std::sqrt(Caloresolution*Caloresolution + muonHCALError + muonECALError);

    // Make PF candidates with the remaining tracks in the block
    sumpError2 = 0.;
    for ( IT it = associatedTracks.begin(); it != associatedTracks.end(); ++it ) { 
      unsigned iTrack = it->second.first;
      if ( !active[iTrack] ) continue;
      reco::TrackRef trackRef = elements[iTrack].trackRef();
      double trackMomentum = trackRef->p();
      double Dp = trackRef->qoverpError()*trackMomentum*trackMomentum;
      unsigned tmpi = reconstructTrack( elements[iTrack] );
      (*pfCandidates_)[tmpi].addElementInBlock( blockref, iTrack );
      (*pfCandidates_)[tmpi].addElementInBlock( blockref, iHcal );
      std::pair<II,II> myEcals = associatedEcals.equal_range(iTrack);
      for (II ii=myEcals.first; ii!=myEcals.second; ++ii ) { 
        unsigned iEcal = ii->second.second;
        if ( active[iEcal] ) continue;
        (*pfCandidates_)[tmpi].addElementInBlock( blockref, iEcal );
      }
      
      if (useHO_) {
        std::pair<II,II> myHOs = associatedHOs.equal_range(iTrack);
        for (II ii=myHOs.first; ii!=myHOs.second; ++ii ) { 
          unsigned iHO = ii->second.second;
          if ( active[iHO] ) continue;
          (*pfCandidates_)[tmpi].addElementInBlock( blockref, iHO );
        }
      }

      if ( iTrack == corrTrack ) { 
        (*pfCandidates_)[tmpi].rescaleMomentum(corrFact);
        trackMomentum *= corrFact;
      }
      chargedHadronsIndices.push_back( tmpi );
      chargedHadronsInBlock.push_back( iTrack );
      active[iTrack] = false;
      hcalP.push_back(trackMomentum);
      hcalDP.push_back(Dp);
      if (Dp/trackMomentum > maxDPovP) maxDPovP = Dp/trackMomentum;
      sumpError2 += Dp*Dp;
    }

    // The total uncertainty of the difference Calo-Track
    TotalError = sqrt(sumpError2 + Caloresolution*Caloresolution);

    if ( debug_ ) {
      cout<<"\tCompare Calo Energy to total charged momentum "<<endl;
      cout<<"\t\tsum p    = "<<totalChargedMomentum<<" +- "<<sqrt(sumpError2)<<endl;
      cout<<"\t\tsum ecal = "<<totalEcal<<endl;
      cout<<"\t\tsum hcal = "<<totalHcal<<endl;
      cout<<"\t\t => Calo Energy = "<<caloEnergy<<" +- "<<Caloresolution<<endl;
      cout<<"\t\t => Calo Energy- total charged momentum = "
          <<caloEnergy-totalChargedMomentum<<" +- "<<TotalError<<endl;
      cout<<endl;
    }

    /* */

    double nsigma = nSigmaHCAL(totalChargedMomentum,hclusterref->positionREP().Eta());
    //double nsigma = nSigmaHCAL(caloEnergy,hclusterref->positionREP().Eta());
    if ( abs(totalChargedMomentum-caloEnergy)<nsigma*TotalError ) {

      // deposited caloEnergy compatible with total charged momentum
      // if tracking errors are large take weighted average
      
      if(debug_) {
        cout<<"\t\tcase 1: COMPATIBLE "
            <<"|Calo Energy- total charged momentum| = "
            <<abs(caloEnergy-totalChargedMomentum)
            <<" < "<<nsigma<<" x "<<TotalError<<endl;
        if (maxDPovP < 0.1 )
          cout<<"\t\t\tmax DP/P = "<< maxDPovP 
              <<" less than 0.1: do nothing "<<endl;
        else 
          cout<<"\t\t\tmax DP/P = "<< maxDPovP 
              <<" >  0.1: take weighted averages "<<endl;
      }
      
      // if max DP/P < 10%  do nothing
      if (maxDPovP > 0.1) {
      
        // for each track associated to hcal
        //      int nrows = tkIs.size();
        int nrows = chargedHadronsIndices.size();
        TMatrixTSym<double> a (nrows);
        TVectorD b (nrows);
        TVectorD check(nrows);
        double sigma2E = Caloresolution*Caloresolution;
        for(int i=0; i<nrows; i++) {
          double sigma2i = hcalDP[i]*hcalDP[i];
          if (debug_){
            cout<<"\t\t\ttrack associated to hcal "<<i 
                <<" P = "<<hcalP[i]<<" +- "
                <<hcalDP[i]<<endl;
          }
          a(i,i) = 1./sigma2i + 1./sigma2E;
          b(i) = hcalP[i]/sigma2i + caloEnergy/sigma2E;
          for(int j=0; j<nrows; j++) {
            if (i==j) continue;
            a(i,j) = 1./sigma2E;
          } // end loop on j
        } // end loop on i
        
        // solve ax = b
        //if (debug_){// debug1
        //cout<<" matrix a(i,j) "<<endl;
        //a.Print();
        //cout<<" vector b(i) "<<endl;
        //b.Print();
        //} // debug1
        TDecompChol decomp(a);
        bool ok = false;
        TVectorD x = decomp.Solve(b, ok);
        // for each track create a PFCandidate track
        // with a momentum rescaled to weighted average
        if (ok) {
          for (int i=0; i<nrows; i++){
            //      unsigned iTrack = trackInfos[i].index;
            unsigned ich = chargedHadronsIndices[i];
            double rescaleFactor =  x(i)/hcalP[i];
            (*pfCandidates_)[ich].rescaleMomentum( rescaleFactor );

            if(debug_){
              cout<<"\t\t\told p "<<hcalP[i]
                  <<" new p "<<x(i) 
                  <<" rescale "<<rescaleFactor<<endl;
            }
          }
        }
        else {
          cerr<<"not ok!"<<endl;
          assert(0);
        }
      }
    }

    else if( caloEnergy > totalChargedMomentum ) {
      
      //case 2: caloEnergy > totalChargedMomentum + nsigma*TotalError
      //there is an excess of energy in the calos
      //create a neutral hadron or a photon

      /*        
      //If it's isolated don't create neutrals since the energy deposit is always coming from a showering muon
      bool thisIsAnIsolatedMuon = false;
      for(IE ie = sortedTracks.begin(); ie != sortedTracks.end(); ++ie ) {
        unsigned iTrack = ie->second;
        if(PFMuonAlgo::isIsolatedMuon(elements[iTrack].muonRef())) thisIsAnIsolatedMuon = true; 
      }
      
      if(thisIsAnIsolatedMuon){
        if(debug_)cout<<" Not looking for neutral b/c this is an isolated muon "<<endl;
        break;
      }
      */
      double eNeutralHadron = caloEnergy - totalChargedMomentum;
      double ePhoton = (caloEnergy - totalChargedMomentum) / slopeEcal;

      if(debug_) {
        if(!sortedTracks.empty() ){
          cout<<"\t\tcase 2: NEUTRAL DETECTION "
              <<caloEnergy<<" > "<<nsigma<<"x"<<TotalError
              <<" + "<<totalChargedMomentum<<endl;
          cout<<"\t\tneutral activity detected: "<<endl
              <<"\t\t\t           photon = "<<ePhoton<<endl
              <<"\t\t\tor neutral hadron = "<<eNeutralHadron<<endl;
            
          cout<<"\t\tphoton or hadron ?"<<endl;}
          
        if(sortedTracks.empty() )
          cout<<"\t\tno track -> hadron "<<endl;
        else 
          cout<<"\t\t"<<sortedTracks.size()
              <<"tracks -> check if the excess is photonic or hadronic"<<endl;
      }
        

      double ratioMax = 0.;
      reco::PFClusterRef maxEcalRef;
      unsigned maxiEcal= 9999;       

      // for each track associated to hcal: iterator IE ie :
        
      for(IE ie = sortedTracks.begin(); ie != sortedTracks.end(); ++ie ) {
          
        unsigned iTrack = ie->second;
          
        PFBlockElement::Type type = elements[iTrack].type();
        assert( type == reco::PFBlockElement::TRACK );
          
        reco::TrackRef trackRef = elements[iTrack].trackRef();
        assert( !trackRef.isNull() ); 
          
        II iae = associatedEcals.find(iTrack);
          
        if( iae == associatedEcals.end() ) continue;
          
        // double distECAL = iae->second.first;
        unsigned iEcal = iae->second.second;
          
        PFBlockElement::Type typeEcal = elements[iEcal].type();
        assert(typeEcal == reco::PFBlockElement::ECAL);
          
        reco::PFClusterRef clusterRef = elements[iEcal].clusterRef();
        assert( !clusterRef.isNull() );
          
        double pTrack = trackRef->p();
        double eECAL = clusterRef->energy();
        double eECALOverpTrack = eECAL / pTrack;

        if ( eECALOverpTrack > ratioMax ) { 
          ratioMax = eECALOverpTrack;
          maxEcalRef = clusterRef;
          maxiEcal = iEcal;
        }          
        
      }  // end loop on tracks associated to hcal: iterator IE ie
    
      std::vector<reco::PFClusterRef> pivotalClusterRef;
      std::vector<unsigned> iPivotal;
      std::vector<double> particleEnergy, ecalEnergy, hcalEnergy, rawecalEnergy, rawhcalEnergy;
      std::vector<math::XYZVector> particleDirection;

      // If the excess is smaller than the ecal energy, assign the whole 
      // excess to a photon
      if ( ePhoton < totalEcal || eNeutralHadron-calibEcal < 1E-10 ) { 
        
        if ( !maxEcalRef.isNull() ) { 
          particleEnergy.push_back(ePhoton);
          particleDirection.push_back(photonAtECAL);
          ecalEnergy.push_back(ePhoton);
          hcalEnergy.push_back(0.);
          rawecalEnergy.push_back(totalEcal);
          rawhcalEnergy.push_back(totalHcal);
          pivotalClusterRef.push_back(maxEcalRef);
          iPivotal.push_back(maxiEcal);
          // So the merged photon energy is,
          mergedPhotonEnergy  = ePhoton;
        }
        
      } else { 
        
        // Otherwise assign the whole ECAL energy to the photon
        if ( !maxEcalRef.isNull() ) { 
          pivotalClusterRef.push_back(maxEcalRef);
          particleEnergy.push_back(totalEcal);
          particleDirection.push_back(photonAtECAL);
          ecalEnergy.push_back(totalEcal);
          hcalEnergy.push_back(0.);
          rawecalEnergy.push_back(totalEcal);
          rawhcalEnergy.push_back(totalHcal);
          iPivotal.push_back(maxiEcal);
          // So the merged photon energy is,
          mergedPhotonEnergy  = totalEcal;
        }
        
        // ... and assign the remaining excess to a neutral hadron
        mergedNeutralHadronEnergy = eNeutralHadron-calibEcal;
        if ( mergedNeutralHadronEnergy > 1.0 ) { 
          pivotalClusterRef.push_back(hclusterref);
          particleEnergy.push_back(mergedNeutralHadronEnergy);
          particleDirection.push_back(hadronAtECAL);
          ecalEnergy.push_back(0.);
          hcalEnergy.push_back(mergedNeutralHadronEnergy);
          rawecalEnergy.push_back(totalEcal);
          rawhcalEnergy.push_back(totalHcal);
          iPivotal.push_back(iHcal);    
        }
        
      }
        

      // reconstructing a merged neutral
      // the type of PFCandidate is known from the 
      // reference to the pivotal cluster. 
      
      for ( unsigned iPivot=0; iPivot<iPivotal.size(); ++iPivot ) { 
        
        if ( particleEnergy[iPivot] < 0. ) 
          std::cout << "ALARM = Negative energy ! " 
                    << particleEnergy[iPivot] << std::endl;
        
        bool useDirection = true;
        unsigned tmpi = reconstructCluster( *pivotalClusterRef[iPivot], 
                                            particleEnergy[iPivot], 
                                            useDirection,
                                            particleDirection[iPivot].X(),
                                            particleDirection[iPivot].Y(),
                                            particleDirection[iPivot].Z()); 

      
        (*pfCandidates_)[tmpi].setEcalEnergy( rawecalEnergy[iPivot],ecalEnergy[iPivot] );
        if ( !useHO_ ) { 
          (*pfCandidates_)[tmpi].setHcalEnergy( rawhcalEnergy[iPivot],hcalEnergy[iPivot] );
          (*pfCandidates_)[tmpi].setHoEnergy(0., 0.);
        } else { 
          (*pfCandidates_)[tmpi].setHcalEnergy( rawhcalEnergy[iPivot]-totalHO,hcalEnergy[iPivot]*(1.-totalHO/rawhcalEnergy[iPivot]));
          (*pfCandidates_)[tmpi].setHoEnergy(totalHO, totalHO * hcalEnergy[iPivot]/rawhcalEnergy[iPivot]);
        } 
        (*pfCandidates_)[tmpi].setPs1Energy( 0. );
        (*pfCandidates_)[tmpi].setPs2Energy( 0. );
        (*pfCandidates_)[tmpi].set_mva_nothing_gamma( -1. );
        //       (*pfCandidates_)[tmpi].addElement(&elements[iPivotal]);
        // (*pfCandidates_)[tmpi].addElementInBlock(blockref, iPivotal[iPivot]);
        (*pfCandidates_)[tmpi].addElementInBlock( blockref, iHcal );
        for ( unsigned ich=0; ich<chargedHadronsInBlock.size(); ++ich) { 
          unsigned iTrack = chargedHadronsInBlock[ich];
          (*pfCandidates_)[tmpi].addElementInBlock( blockref, iTrack );
          // Assign the position of the track at the ECAL entrance
          const math::XYZPointF& chargedPosition = 
            dynamic_cast<const reco::PFBlockElementTrack*>(&elements[iTrack])->positionAtECALEntrance();
          (*pfCandidates_)[tmpi].setPositionAtECALEntrance(chargedPosition);

          std::pair<II,II> myEcals = associatedEcals.equal_range(iTrack);
          for (II ii=myEcals.first; ii!=myEcals.second; ++ii ) { 
            unsigned iEcal = ii->second.second;
            if ( active[iEcal] ) continue;
            (*pfCandidates_)[tmpi].addElementInBlock( blockref, iEcal );
          }
        }

      }

    } // excess of energy
  
  
    // will now share the hcal energy between the various charged hadron 
    // candidates, taking into account the potential neutral hadrons
    
    double totalHcalEnergyCalibrated = calibHcal;
    double totalEcalEnergyCalibrated = calibEcal;
    //JB: The question is: we've resolved the merged photons cleanly, but how should
    //the remaining hadrons be assigned the remaining ecal energy?
    //*Temporary solution*: follow HCAL example with fractions...    
 
   /*
    if(totalEcal>0) {
      // removing ecal energy from abc calibration
      totalHcalEnergyCalibrated  -= calibEcal;
      // totalHcalEnergyCalibrated -= calibration_->paramECALplusHCAL_slopeECAL() * totalEcal;
    }
    */
    // else caloEnergy = hcal only calibrated energy -> ok.

    // remove the energy of the potential neutral hadron
    totalHcalEnergyCalibrated -= mergedNeutralHadronEnergy;
    // similarly for the merged photons
    totalEcalEnergyCalibrated -= mergedPhotonEnergy;
    // share between the charged hadrons

    //COLIN can compute this before
    // not exactly equal to sum p, this is sum E
    double chargedHadronsTotalEnergy = 0;
    for( unsigned ich=0; ich<chargedHadronsIndices.size(); ++ich ) {
      unsigned index = chargedHadronsIndices[ich];
      reco::PFCandidate& chargedHadron = (*pfCandidates_)[index];
      chargedHadronsTotalEnergy += chargedHadron.energy();
    }

    for( unsigned ich=0; ich<chargedHadronsIndices.size(); ++ich ) {
      unsigned index = chargedHadronsIndices[ich];
      reco::PFCandidate& chargedHadron = (*pfCandidates_)[index];
      float fraction = chargedHadron.energy()/chargedHadronsTotalEnergy;

      if ( !useHO_ ) { 
        chargedHadron.setHcalEnergy(  fraction * totalHcal, fraction * totalHcalEnergyCalibrated );          
        chargedHadron.setHoEnergy(  0., 0. ); 
      } else { 
        chargedHadron.setHcalEnergy(  fraction * (totalHcal-totalHO), fraction * totalHcalEnergyCalibrated * (1.-totalHO/totalHcal) );          
        chargedHadron.setHoEnergy( fraction * totalHO, fraction * totalHO * totalHcalEnergyCalibrated / totalHcal ); 
      }
      //JB: fixing up (previously omitted) setting of ECAL energy gouzevit
      chargedHadron.setEcalEnergy( fraction * totalEcal, fraction * totalEcalEnergyCalibrated );
    }

    // Finally treat unused ecal satellites as photons.
    for ( IS is = ecalSatellites.begin(); is != ecalSatellites.end(); ++is ) { 
      
      // Ignore satellites already taken
      unsigned iEcal = is->second.first;
      if ( !active[iEcal] ) continue;

      // Sanity checks again (well not useful, this time!)
      PFBlockElement::Type type = elements[ iEcal ].type();
      assert( type == PFBlockElement::ECAL );
      PFClusterRef eclusterref = elements[iEcal].clusterRef();
      assert(!eclusterref.isNull() );

      // Lock the cluster
      active[iEcal] = false;
      
      // Find the associated tracks
      std::multimap<double, unsigned> assTracks;
      block.associatedElements( iEcal,  linkData,
                                assTracks,
                                reco::PFBlockElement::TRACK,
                                reco::PFBlock::LINKTEST_ALL );

      // Create a photon
      unsigned tmpi = reconstructCluster( *eclusterref, sqrt(is->second.second.Mag2()) ); 
      (*pfCandidates_)[tmpi].setEcalEnergy( eclusterref->energy(),sqrt(is->second.second.Mag2()) );
      (*pfCandidates_)[tmpi].setHcalEnergy( 0., 0. );
      (*pfCandidates_)[tmpi].setHoEnergy( 0., 0. );
      (*pfCandidates_)[tmpi].setPs1Energy( associatedPSs[iEcal].first );
      (*pfCandidates_)[tmpi].setPs2Energy( associatedPSs[iEcal].second );
      (*pfCandidates_)[tmpi].addElementInBlock( blockref, iEcal );
      (*pfCandidates_)[tmpi].addElementInBlock( blockref, sortedTracks.begin()->second) ;
    }


  }

   // end loop on hcal element iHcal= hcalIs[i] 


  // Processing the remaining HCAL clusters
  if(debug_) { 
    cout<<endl;
    if(debug_) 
      cout<<endl
          <<"---- loop remaining hcal ------- "<<endl;
  }
  

  //COLINFEB16 
  // now dealing with the HCAL elements that are not linked to any track 
  for(unsigned ihcluster=0; ihcluster<hcalIs.size(); ihcluster++) {
    unsigned iHcal = hcalIs[ihcluster];

    // Keep ECAL and HO elements for reference in the PFCandidate
    std::vector<unsigned> ecalRefs;
    std::vector<unsigned> hoRefs;
    
    if(debug_)
      cout<<endl<<elements[iHcal]<<" ";
    
    
    if( !active[iHcal] ) {
      if(debug_) 
        cout<<"not active"<<endl;         
      continue;
    }
    
    // Find the ECAL elements linked to it
    std::multimap<double, unsigned> ecalElems;
    block.associatedElements( iHcal,  linkData,
                              ecalElems ,
                              reco::PFBlockElement::ECAL,
                              reco::PFBlock::LINKTEST_ALL );

    // Loop on these ECAL elements
    float totalEcal = 0.;
    float ecalMax = 0.;
    reco::PFClusterRef eClusterRef;
    for(IE ie = ecalElems.begin(); ie != ecalElems.end(); ++ie ) {
      
      unsigned iEcal = ie->second;
      double dist = ie->first;
      PFBlockElement::Type type = elements[iEcal].type();
      assert( type == PFBlockElement::ECAL );
      
      // Check if already used
      if( !active[iEcal] ) continue;

      // Check the distance (one HCALPlusECAL tower, roughly)
      // if ( dist > 0.15 ) continue;

      //COLINFEB16 
      // what could be done is to
      // - link by rechit.
      // - take in the neutral hadron all the ECAL clusters 
      // which are within the same CaloTower, according to the distance, 
      // except the ones which are closer to another HCAL cluster.
      // - all the other ECAL linked to this HCAL are photons. 
      // 
      // about the closest HCAL cluster. 
      // it could maybe be easier to loop on the ECAL clusters first
      // to cut the links to all HCAL clusters except the closest, as is 
      // done in the first track loop. But maybe not!
      // or add an helper function to the PFAlgo class to ask 
      // if a given element is the closest of a given type to another one?
    
      // Check if not closer from another free HCAL 
      std::multimap<double, unsigned> hcalElems;
      block.associatedElements( iEcal,  linkData,
                                hcalElems ,
                                reco::PFBlockElement::HCAL,
                                reco::PFBlock::LINKTEST_ALL );

      bool isClosest = true;
      for(IE ih = hcalElems.begin(); ih != hcalElems.end(); ++ih ) {

        unsigned jHcal = ih->second;
        double distH = ih->first;
        
        if ( !active[jHcal] ) continue;

        if ( distH < dist ) { 
          isClosest = false;
          break;
        }

      }

      if (!isClosest) continue;


      if(debug_) {
        cout<<"\telement "<<elements[iEcal]<<" linked with dist "<< dist<<endl;
        cout<<"Added to HCAL cluster to form a neutral hadron"<<endl;
      }
      
      reco::PFClusterRef eclusterRef = elements[iEcal].clusterRef();
      assert( !eclusterRef.isNull() );

      // Check the presence of ps clusters in the vicinity
      vector<double> ps1Ene(1,static_cast<double>(0.));
      associatePSClusters(iEcal, reco::PFBlockElement::PS1, block, elements, linkData, active, ps1Ene);
      vector<double> ps2Ene(1,static_cast<double>(0.));
      associatePSClusters(iEcal, reco::PFBlockElement::PS2, block, elements, linkData, active, ps2Ene);
      bool crackCorrection = false;
      double ecalEnergy = calibration_->energyEm(*eclusterRef,ps1Ene,ps2Ene,crackCorrection);
      
      //std::cout << "EcalEnergy, ps1, ps2 = " << ecalEnergy 
      //          << ", " << ps1Ene[0] << ", " << ps2Ene[0] << std::endl;
      totalEcal += ecalEnergy;
      if ( ecalEnergy > ecalMax ) { 
        ecalMax = ecalEnergy;
        eClusterRef = eclusterRef;
      }
      
      ecalRefs.push_back(iEcal);
      active[iEcal] = false;


    }// End loop ECAL
    
    // Now find the HO clusters linked to the HCAL cluster
    double totalHO = 0.;
    double hoMax = 0.;
    //unsigned jHO = 0;
    if (useHO_) {
      std::multimap<double, unsigned> hoElems;
      block.associatedElements( iHcal,  linkData,
                                hoElems ,
                                reco::PFBlockElement::HO,
                                reco::PFBlock::LINKTEST_ALL );
      
      // Loop on these HO elements
      //      double totalHO = 0.;
      //      double hoMax = 0.;
      //      unsigned jHO = 0;
      reco::PFClusterRef hoClusterRef;
      for(IE ie = hoElems.begin(); ie != hoElems.end(); ++ie ) {
        
        unsigned iHO = ie->second;
        double dist = ie->first;
        PFBlockElement::Type type = elements[iHO].type();
        assert( type == PFBlockElement::HO );
        
        // Check if already used
        if( !active[iHO] ) continue;
        
        // Check the distance (one HCALPlusHO tower, roughly)
        // if ( dist > 0.15 ) continue;
        
        // Check if not closer from another free HCAL 
        std::multimap<double, unsigned> hcalElems;
        block.associatedElements( iHO,  linkData,
                                  hcalElems ,
                                  reco::PFBlockElement::HCAL,
                                  reco::PFBlock::LINKTEST_ALL );
        
        bool isClosest = true;
        for(IE ih = hcalElems.begin(); ih != hcalElems.end(); ++ih ) {
          
          unsigned jHcal = ih->second;
          double distH = ih->first;
          
          if ( !active[jHcal] ) continue;
          
          if ( distH < dist ) { 
            isClosest = false;
            break;
          }
          
        }
        
        if (!isClosest) continue;
        
        if(debug_ && useHO_) {
          cout<<"\telement "<<elements[iHO]<<" linked with dist "<< dist<<endl;
          cout<<"Added to HCAL cluster to form a neutral hadron"<<endl;
        }
        
        reco::PFClusterRef hoclusterRef = elements[iHO].clusterRef();
        assert( !hoclusterRef.isNull() );
        
        double hoEnergy = hoclusterRef->energy(); // calibration_->energyEm(*hoclusterRef,ps1Ene,ps2Ene,crackCorrection);
        
        totalHO += hoEnergy;
        if ( hoEnergy > hoMax ) { 
          hoMax = hoEnergy;
          hoClusterRef = hoclusterRef;
          //jHO = iHO;
        }
        
        hoRefs.push_back(iHO);
        active[iHO] = false;
        
      }// End loop HO
    }

    PFClusterRef hclusterRef 
      = elements[iHcal].clusterRef();
    assert( !hclusterRef.isNull() );
    
    // HCAL energy
    double totalHcal =  hclusterRef->energy();
    // Include the HO energy 
    if ( useHO_ ) totalHcal += totalHO;

    // std::cout << "Hcal Energy,eta : " << totalHcal 
    //          << ", " << hclusterRef->positionREP().Eta()
    //          << std::endl;
    // Calibration
    //double caloEnergy = totalHcal;
    // double slopeEcal = 1.0;
    double calibEcal = totalEcal > 0. ? totalEcal : 0.;
    double calibHcal = std::max(0.,totalHcal);
    if  (  hclusterRef->layer() == PFLayer::HF_HAD  ||
           hclusterRef->layer() == PFLayer::HF_EM ) { 
      //caloEnergy = totalHcal/0.7;
      calibEcal = totalEcal;
    } else { 
      calibration_->energyEmHad(-1.,calibEcal,calibHcal,
                                hclusterRef->positionREP().Eta(),
                                hclusterRef->positionREP().Phi());      
      //caloEnergy = calibEcal+calibHcal;
    }

    // std::cout << "CalibEcal,HCal = " << calibEcal << ", " << calibHcal << std::endl;
    // std::cout << "-------------------------------------------------------------------" << std::endl;
    // ECAL energy : calibration

    // double particleEnergy = caloEnergy;
    // double particleEnergy = totalEcal + calibHcal;
    // particleEnergy /= (1.-0.724/sqrt(particleEnergy)-0.0226/particleEnergy);

    unsigned tmpi = reconstructCluster( *hclusterRef, 
                                        calibEcal+calibHcal ); 

    
    (*pfCandidates_)[tmpi].setEcalEnergy( totalEcal, calibEcal );
    if ( !useHO_ ) { 
      (*pfCandidates_)[tmpi].setHcalEnergy( totalHcal, calibHcal );
      (*pfCandidates_)[tmpi].setHoEnergy(0.,0.);
    } else { 
      (*pfCandidates_)[tmpi].setHcalEnergy( totalHcal-totalHO, calibHcal*(1.-totalHO/totalHcal));
      (*pfCandidates_)[tmpi].setHoEnergy(totalHO,totalHO*calibHcal/totalHcal);
    }
    (*pfCandidates_)[tmpi].setPs1Energy( 0. );
    (*pfCandidates_)[tmpi].setPs2Energy( 0. );
    (*pfCandidates_)[tmpi].addElementInBlock( blockref, iHcal );
    for (unsigned iec=0; iec<ecalRefs.size(); ++iec) 
      (*pfCandidates_)[tmpi].addElementInBlock( blockref, ecalRefs[iec] );
    for (unsigned iho=0; iho<hoRefs.size(); ++iho) 
      (*pfCandidates_)[tmpi].addElementInBlock( blockref, hoRefs[iho] );
      
  }//loop hcal elements




  if(debug_) { 
    cout<<endl;
    if(debug_) cout<<endl<<"---- loop ecal------- "<<endl;
  }
  
  // for each ecal element iEcal = ecalIs[i] in turn:

  for(unsigned i=0; i<ecalIs.size(); i++) {
    unsigned iEcal = ecalIs[i];
    
    if(debug_) 
      cout<<endl<<elements[iEcal]<<" ";
    
    if( ! active[iEcal] ) {
      if(debug_) 
        cout<<"not active"<<endl;         
      continue;
    }
    
    PFBlockElement::Type type = elements[ iEcal ].type();
    assert( type == PFBlockElement::ECAL );

    PFClusterRef clusterref = elements[iEcal].clusterRef();
    assert(!clusterref.isNull() );

    active[iEcal] = false;

    // Check the presence of ps clusters in the vicinity
    vector<double> ps1Ene(1,static_cast<double>(0.));
    associatePSClusters(iEcal, reco::PFBlockElement::PS1, block, elements, linkData, active, ps1Ene);
    vector<double> ps2Ene(1,static_cast<double>(0.));
    associatePSClusters(iEcal, reco::PFBlockElement::PS2, block, elements, linkData, active, ps2Ene);
    bool crackCorrection = false;
    float ecalEnergy = calibration_->energyEm(*clusterref,ps1Ene,ps2Ene,crackCorrection);
    // float ecalEnergy = calibration_->energyEm( clusterref->energy() );
    double particleEnergy = ecalEnergy;
    
    unsigned tmpi = reconstructCluster( *clusterref, 
                                        particleEnergy );
 
    (*pfCandidates_)[tmpi].setEcalEnergy( clusterref->energy(),ecalEnergy );
    (*pfCandidates_)[tmpi].setHcalEnergy( 0., 0. );
    (*pfCandidates_)[tmpi].setHoEnergy( 0., 0. );
    (*pfCandidates_)[tmpi].setPs1Energy( 0. );
    (*pfCandidates_)[tmpi].setPs2Energy( 0. );
    (*pfCandidates_)[tmpi].addElementInBlock( blockref, iEcal );
    

  }  // end loop on ecal elements iEcal = ecalIs[i]

}  // end processBlock



unsigned PFAlgo::reconstructTrack( const reco::PFBlockElement& elt ) {

  const reco::PFBlockElementTrack* eltTrack 
    = dynamic_cast<const reco::PFBlockElementTrack*>(&elt);

  reco::TrackRef trackRef = eltTrack->trackRef();
  const reco::Track& track = *trackRef;

  reco::MuonRef muonRef = eltTrack->muonRef();

  int charge = track.charge()>0 ? 1 : -1;

  // Assign the pion mass to all charged particles
  double px = track.px();
  double py = track.py();
  double pz = track.pz();
  double energy = sqrt(track.p()*track.p() + 0.13957*0.13957);


  // Except if it is a muon, of course ! 
  bool thisIsAMuon = PFMuonAlgo::isMuon(elt);
  bool thisIsAnIsolatedMuon = PFMuonAlgo::isIsolatedMuon(elt);
  bool thisIsAGlobalTightMuon = PFMuonAlgo::isGlobalTightMuon(elt);
  bool thisIsATrackerTightMuon = PFMuonAlgo::isTrackerTightMuon(elt);

  // Or from nuclear inetraction then use the refitted momentum
  bool isFromDisp = isFromSecInt(elt, "secondary");
  bool isToDisp  = isFromSecInt(elt, "primary");
  //isFromNucl = false;
  bool globalFitUsed = false;

  if ( thisIsAMuon ) { 
    
    //By default take muon kinematics directly from the muon object (usually determined by the track) 
    px = muonRef->px();
    py = muonRef->py();
    pz = muonRef->pz();
    energy = sqrt(muonRef->p()*muonRef->p() + 0.1057*0.1057);

    reco::TrackBase::TrackQuality trackQualityHighPurity = TrackBase::qualityByName("highPurity");
    if(debug_)if(!trackRef->quality(trackQualityHighPurity))cout<<" Low Purity Track "<<endl;

    if(muonRef->isGlobalMuon() && !usePFMuonMomAssign_){
      
      reco::TrackRef bestMuTrack =  muonRef->muonBestTrack();
      if(!useBestMuonTrack_)
        bestMuTrack = muonRef->combinedMuon(); 

      // take the global fit instead under special circumstances
      bool useGlobalFit = false;
      
      if(thisIsAnIsolatedMuon && (!muonRef->isTrackerMuon() || (muonRef->pt() > bestMuTrack->pt() && track.ptError() > 5.0*bestMuTrack->ptError()))) useGlobalFit = true;
      else if(!trackRef->quality(trackQualityHighPurity)) useGlobalFit = true;
      else if(muonRef->pt() > bestMuTrack->pt() &&
              (track.hitPattern().numberOfValidTrackerHits() < 8 || track.hitPattern().numberOfValidPixelHits() == 0 ) &&
              track.ptError() > 5.0*bestMuTrack->ptError()) useGlobalFit = true;

      if(useGlobalFit){
        px = bestMuTrack->px();
        py = bestMuTrack->py();
        pz = bestMuTrack->pz();
        energy = sqrt(bestMuTrack->p()*bestMuTrack->p() + 0.1057*0.1057);   
        globalFitUsed = true;
      }
    }      
    else if(usePFMuonMomAssign_){
      // If this option is set we take more liberties choosing the muon kinematics (not advised by the muon POG)
      if(thisIsAGlobalTightMuon)
        {
          // If the global muon above 10 GeV and is a tracker muon take the global pT     
          if(muonRef->isTrackerMuon()){
            if(sqrt(px*px+py*py) > 10){

              reco::TrackRef bestMuTrack =  muonRef->muonBestTrack();
              if(!useBestMuonTrack_)
                bestMuTrack = muonRef->combinedMuon(); 
              
              px = bestMuTrack->px();
              py = bestMuTrack->py();
              pz = bestMuTrack->pz();
              energy = sqrt(bestMuTrack->p()*bestMuTrack->p() + 0.1057*0.1057);   
              globalFitUsed = true;
            }
          }   // If it's not a tracker muon, choose between the global pT and the STA pT
          else{ 
            
            reco::TrackRef bestMuTrack = muonRef->combinedMuon()->normalizedChi2() < muonRef->standAloneMuon()->normalizedChi2() ?
              muonRef->muonBestTrack() : 
              muonRef->standAloneMuon() ;
            
            if(!useBestMuonTrack_)
              bestMuTrack =
                muonRef->combinedMuon()->normalizedChi2() < muonRef->standAloneMuon()->normalizedChi2() ?
                muonRef->combinedMuon() : 
                muonRef->standAloneMuon() ;

            px = bestMuTrack->px();
            py = bestMuTrack->py();
            pz = bestMuTrack->pz();
            energy = sqrt(bestMuTrack->p()*bestMuTrack->p() + 0.1057*0.1057);   
            globalFitUsed = true;
          }
        } // close else if(thisIsAGlobalTightMuon)
    } // close (usePFPFMuonMomAssign_)      
  }// close if(thisIsAMuon)
  else if (isFromDisp) {
    double Dpt = trackRef->ptError();
    double dptRel = Dpt/trackRef->pt()*100;
    //If the track is ill measured it is better to not refit it, since the track information probably would not be used.
    //In the PFAlgo we use the trackref information. If the track error is too big the refitted information might be very different
    // from the not refitted one.
    if (dptRel < dptRel_DispVtx_){

      if (debug_) 
        cout << "Not refitted px = " << px << " py = " << py << " pz = " << pz << " energy = " << energy << endl; 
      //reco::TrackRef trackRef = eltTrack->trackRef();
      reco::PFDisplacedVertexRef vRef = eltTrack->displacedVertexRef(reco::PFBlockElement::T_FROM_DISP)->displacedVertexRef();
      reco::Track trackRefit = vRef->refittedTrack(trackRef);
      px = trackRefit.px();
      py = trackRefit.py();
      pz = trackRefit.pz();
      energy = sqrt(trackRefit.p()*trackRefit.p() + 0.13957*0.13957);
      if (debug_) 
        cout << "Refitted px = " << px << " py = " << py << " pz = " << pz << " energy = " << energy << endl; 
    
    }
  }
  
  if ((isFromDisp || isToDisp) && debug_) cout << "Final px = " << px << " py = " << py << " pz = " << pz << " energy = " << energy << endl; 

  // Create a PF Candidate
  math::XYZTLorentzVector momentum(px,py,pz,energy);
  reco::PFCandidate::ParticleType particleType 
    = reco::PFCandidate::h;

  // Add it to the stack
  pfCandidates_->push_back( PFCandidate( charge, 
                                         momentum,
                                         particleType ) );

  // displaced vertices 
  if( isFromDisp ) {
    pfCandidates_->back().setFlag( reco::PFCandidate::T_FROM_DISP, true);
    pfCandidates_->back().setDisplacedVertexRef( eltTrack->displacedVertexRef(reco::PFBlockElement::T_FROM_DISP)->displacedVertexRef(), reco::PFCandidate::T_FROM_DISP);
  }
  
  // do not label as primary a track which would be recognised as a muon. A muon cannot produce NI. It is with high probability a fake
  if( isToDisp && !thisIsAMuon ) {
    pfCandidates_->back().setFlag( reco::PFCandidate::T_TO_DISP, true);
    pfCandidates_->back().setDisplacedVertexRef( eltTrack->displacedVertexRef(reco::PFBlockElement::T_TO_DISP)->displacedVertexRef(), reco::PFCandidate::T_TO_DISP);
  }



  if ( thisIsAMuon && globalFitUsed ) 
    pfCandidates_->back().setVertexSource(  PFCandidate::kComMuonVertex );
  else
    pfCandidates_->back().setVertexSource( PFCandidate::kTrkVertex );

  pfCandidates_->back().setTrackRef( trackRef );
  pfCandidates_->back().setPositionAtECALEntrance( eltTrack->positionAtECALEntrance());



  // setting the muon ref if there is
  if (muonRef.isNonnull()) {
    pfCandidates_->back().setMuonRef( muonRef );
    // setting the muon particle type if it is a global muon
    if ( thisIsAMuon) {
      particleType = reco::PFCandidate::mu;
      pfCandidates_->back().setParticleType( particleType );
      if (debug_) {
        if(thisIsAGlobalTightMuon) cout << "PFAlgo: particle type set to muon (global, tight), pT = " <<muonRef->pt()<< endl; 
        else if(thisIsATrackerTightMuon) cout << "PFAlgo: particle type set to muon (tracker, tight), pT = " <<muonRef->pt()<< endl;
        else if(thisIsAnIsolatedMuon) cout << "PFAlgo: particle type set to muon (isolated), pT = " <<muonRef->pt()<< endl; 
        else cout<<" problem with muon assignment "<<endl;
        PFMuonAlgo::printMuonProperties( muonRef );
      }
    }
  }  

  // conversion...

  if(debug_) 
    cout<<"** candidate: "<<pfCandidates_->back()<<endl; 
  
  // returns index to the newly created PFCandidate
  return pfCandidates_->size()-1;
}





unsigned 
PFAlgo::reconstructCluster(const reco::PFCluster& cluster,
                           double particleEnergy, 
                           bool useDirection, 
                           double particleX,
                           double particleY, 
                           double particleZ) {
  
  reco::PFCandidate::ParticleType particleType = reco::PFCandidate::X;

  // need to convert the math::XYZPoint data member of the PFCluster class=
  // to a displacement vector: 

  // Transform particleX,Y,Z to a position at ECAL/HCAL entrance
  double factor = 1.;
  if ( useDirection ) { 
    switch( cluster.layer() ) {
    case PFLayer::ECAL_BARREL:
    case PFLayer::HCAL_BARREL1:
      factor = std::sqrt(cluster.position().Perp2()/(particleX*particleX+particleY*particleY));
      break;
    case PFLayer::ECAL_ENDCAP:
    case PFLayer::HCAL_ENDCAP:
    case PFLayer::HF_HAD:
    case PFLayer::HF_EM:
      factor = cluster.position().Z()/particleZ;
      break;
    default:
      assert(0);
    }
  }
  //MIKE First of all let's check if we have vertex.
  math::XYZPoint vertexPos; 
  if(useVertices_)
    vertexPos = math::XYZPoint(primaryVertex_.x(),primaryVertex_.y(),primaryVertex_.z());
  else
    vertexPos = math::XYZPoint(0.0,0.0,0.0);


  math::XYZVector clusterPos( cluster.position().X()-vertexPos.X(), 
                              cluster.position().Y()-vertexPos.Y(),
                              cluster.position().Z()-vertexPos.Z());
  math::XYZVector particleDirection ( particleX*factor-vertexPos.X(), 
                                      particleY*factor-vertexPos.Y(), 
                                      particleZ*factor-vertexPos.Z() );

  //math::XYZVector clusterPos( cluster.position().X(), cluster.position().Y(),cluster.position().Z() );
  //math::XYZVector particleDirection ( particleX, particleY, particleZ );

  clusterPos = useDirection ? particleDirection.Unit() : clusterPos.Unit();
  clusterPos *= particleEnergy;

  // clusterPos is now a vector along the cluster direction, 
  // with a magnitude equal to the cluster energy.
  
  double mass = 0;
  ROOT::Math::LorentzVector<ROOT::Math::PxPyPzM4D<double> > 
    momentum( clusterPos.X(), clusterPos.Y(), clusterPos.Z(), mass); 
  // mathcore is a piece of #$%
  math::XYZTLorentzVector  tmp;
  // implicit constructor not allowed
  tmp = momentum;

  // Charge
  int charge = 0;

  // Type
  switch( cluster.layer() ) {
  case PFLayer::ECAL_BARREL:
  case PFLayer::ECAL_ENDCAP:
    particleType = PFCandidate::gamma;
    break;
  case PFLayer::HCAL_BARREL1:
  case PFLayer::HCAL_ENDCAP:    
    particleType = PFCandidate::h0;
    break;
  case PFLayer::HF_HAD:
    particleType = PFCandidate::h_HF;
    break;
  case PFLayer::HF_EM:
    particleType = PFCandidate::egamma_HF;
    break;
  default:
    assert(0);
  }

  // The pf candidate
  pfCandidates_->push_back( PFCandidate( charge, 
                                         tmp, 
                                         particleType ) );

  // The position at ECAL entrance (well: watch out, it is not true
  // for HCAL clusters... to be fixed)
  pfCandidates_->back().
    setPositionAtECALEntrance(math::XYZPointF(cluster.position().X(),
                                              cluster.position().Y(),
                                              cluster.position().Z()));

  //Set the cnadidate Vertex
  pfCandidates_->back().setVertex(vertexPos);  

  if(debug_) 
    cout<<"** candidate: "<<pfCandidates_->back()<<endl; 

  // returns index to the newly created PFCandidate
  return pfCandidates_->size()-1;

}


//GMA need the followign two for HO also

double 
PFAlgo::neutralHadronEnergyResolution(double clusterEnergyHCAL, double eta) const {

  // Add a protection
  if ( clusterEnergyHCAL < 1. ) clusterEnergyHCAL = 1.;

  double resol =  fabs(eta) < 1.48 ? 
    sqrt (1.02*1.02/clusterEnergyHCAL + 0.065*0.065)
    :
    sqrt (1.20*1.20/clusterEnergyHCAL + 0.028*0.028);

  return resol;

}

double
PFAlgo::nSigmaHCAL(double clusterEnergyHCAL, double eta) const {
  double nS = fabs(eta) < 1.48 ? 
    nSigmaHCAL_ * (1. + exp(-clusterEnergyHCAL/100.))     
    :
    nSigmaHCAL_ * (1. + exp(-clusterEnergyHCAL/100.));     
  
  return nS;
}


ostream& operator<<(ostream& out, const PFAlgo& algo) {
  if(!out ) return out;

  out<<"====== Particle Flow Algorithm ======= ";
  out<<endl;
  out<<"nSigmaECAL_     "<<algo.nSigmaECAL_<<endl;
  out<<"nSigmaHCAL_     "<<algo.nSigmaHCAL_<<endl;
  out<<endl;
  out<<(*(algo.calibration_))<<endl;
  out<<endl;
  out<<"reconstructed particles: "<<endl;
   
  const std::auto_ptr< reco::PFCandidateCollection >& 
    candidates = algo.pfCandidates(); 

  if(!candidates.get() ) {
    out<<"candidates already transfered"<<endl;
    return out;
  }
  for(PFCandidateConstIterator ic=algo.pfCandidates_->begin(); 
      ic != algo.pfCandidates_->end(); ic++) {
    out<<(*ic)<<endl;
  }

  return out;
}

reco::PFBlockRef 
PFAlgo::createBlockRef( const reco::PFBlockCollection& blocks, 
                        unsigned bi ) {

  if( blockHandle_.isValid() ) {
    return reco::PFBlockRef(  blockHandle_, bi );
  } 
  else {
    return reco::PFBlockRef(  &blocks, bi );
  }
}

void 
PFAlgo::associatePSClusters(unsigned iEcal,
                            reco::PFBlockElement::Type psElementType,
                            const reco::PFBlock& block,
                            const edm::OwnVector< reco::PFBlockElement >& elements, 
                            const reco::PFBlock::LinkData& linkData,
                            std::vector<bool>& active, 
                            std::vector<double>& psEne) {

  // Find all PS clusters with type psElement associated to ECAL cluster iEcal, 
  // within all PFBlockElement "elements" of a given PFBlock "block"
  // psElement can be reco::PFBlockElement::PS1 or reco::PFBlockElement::PS2
  // Returns a vector of PS cluster energies, and updates the "active" vector.

  // Find all PS clusters linked to the iEcal cluster
  std::multimap<double, unsigned> sortedPS;
  typedef std::multimap<double, unsigned>::iterator IE;
  block.associatedElements( iEcal,  linkData,
                            sortedPS, psElementType,
                            reco::PFBlock::LINKTEST_ALL );

  // Loop over these PS clusters
  double totalPS = 0.;
  for ( IE ips=sortedPS.begin(); ips!=sortedPS.end(); ++ips ) {

    // CLuster index and distance to iEcal
    unsigned iPS = ips->second;
    // double distPS = ips->first;

    // Ignore clusters already in use
    if (!active[iPS]) continue;

    // Check that this cluster is not closer to another ECAL cluster
    std::multimap<double, unsigned> sortedECAL;
    block.associatedElements( iPS,  linkData,
                              sortedECAL,
                              reco::PFBlockElement::ECAL,
                              reco::PFBlock::LINKTEST_ALL );
    unsigned jEcal = sortedECAL.begin()->second;
    if ( jEcal != iEcal ) continue; 
    
    // Update PS energy
    PFBlockElement::Type pstype = elements[ iPS ].type();
    assert( pstype == psElementType );
    PFClusterRef psclusterref = elements[iPS].clusterRef();
    assert(!psclusterref.isNull() );
    totalPS += psclusterref->energy(); 
    psEne[0] += psclusterref->energy();
    active[iPS] = false;
  }
            

}


bool
PFAlgo::isFromSecInt(const reco::PFBlockElement& eTrack, string order) const {

  reco::PFBlockElement::TrackType T_TO_DISP = reco::PFBlockElement::T_TO_DISP;
  reco::PFBlockElement::TrackType T_FROM_DISP = reco::PFBlockElement::T_FROM_DISP;
  //  reco::PFBlockElement::TrackType T_FROM_GAMMACONV = reco::PFBlockElement::T_FROM_GAMMACONV;
  reco::PFBlockElement::TrackType T_FROM_V0 = reco::PFBlockElement::T_FROM_V0;

  bool bPrimary = (order.find("primary") != string::npos);
  bool bSecondary = (order.find("secondary") != string::npos);
  bool bAll = (order.find("all") != string::npos);

  bool isToDisp = usePFNuclearInteractions_ && eTrack.trackType(T_TO_DISP);
  bool isFromDisp = usePFNuclearInteractions_ && eTrack.trackType(T_FROM_DISP);

  if (bPrimary && isToDisp) return true;
  if (bSecondary && isFromDisp ) return true;
  if (bAll && ( isToDisp || isFromDisp ) ) return true;

//   bool isFromConv = usePFConversions_ && eTrack.trackType(T_FROM_GAMMACONV);

//   if ((bAll || bSecondary)&& isFromConv) return true;

  bool isFromDecay = (bAll || bSecondary) && usePFDecays_ && eTrack.trackType(T_FROM_V0);

  return isFromDecay;


}


void  
PFAlgo::setEGElectronCollection(const reco::GsfElectronCollection & egelectrons) {
  if(useEGElectrons_ && pfele_) pfele_->setEGElectronCollection(egelectrons);
}

void 
PFAlgo::postCleaning() { 
  
  // Check if the post HF Cleaning was requested - if not, do nothing
  if ( !postHFCleaning_ ) return;

  //Compute met and met significance (met/sqrt(SumEt))
  double metX = 0.;
  double metY = 0.;
  double sumet = 0;
  std::vector<unsigned int> pfCandidatesToBeRemoved;
  for(unsigned i=0; i<pfCandidates_->size(); i++) {
    const PFCandidate& pfc = (*pfCandidates_)[i];
    metX += pfc.px();
    metY += pfc.py();
    sumet += pfc.pt();
  }    
  double met2 = metX*metX+metY*metY;
  // Select events with large MET significance.
  double significance = std::sqrt(met2/sumet);
  double significanceCor = significance;
  if ( significance > minSignificance_ ) { 
    
    double metXCor = metX;
    double metYCor = metY;
    double sumetCor = sumet;
    double met2Cor = met2;
    double deltaPhi = 3.14159;
    double deltaPhiPt = 100.;
    bool next = true;
    unsigned iCor = 1E9;

    // Find the HF candidate with the largest effect on the MET
    while ( next ) { 

      double metReduc = -1.;
      // Loop on the candidates
      for(unsigned i=0; i<pfCandidates_->size(); ++i) {
        const PFCandidate& pfc = (*pfCandidates_)[i];

        // Check that the pfCandidate is in the HF
        if ( pfc.particleId() != reco::PFCandidate::h_HF && 
             pfc.particleId() != reco::PFCandidate::egamma_HF ) continue;

        // Check if has meaningful pt
        if ( pfc.pt() < minHFCleaningPt_ ) continue;
        
        // Check that it is  not already scheculed to be cleaned
        bool skip = false;
        for(unsigned j=0; j<pfCandidatesToBeRemoved.size(); ++j) {
          if ( i == pfCandidatesToBeRemoved[j] ) skip = true;
          if ( skip ) break;
        }
        if ( skip ) continue;
        
        // Check that the pt and the MET are aligned
        deltaPhi = std::acos((metX*pfc.px()+metY*pfc.py())/(pfc.pt()*std::sqrt(met2)));
        deltaPhiPt = deltaPhi*pfc.pt();
        if ( deltaPhiPt > maxDeltaPhiPt_ ) continue;

        // Now remove the candidate from the MET
        double metXInt = metX - pfc.px();
        double metYInt = metY - pfc.py();
        double sumetInt = sumet - pfc.pt();
        double met2Int = metXInt*metXInt+metYInt*metYInt;
        if ( met2Int < met2Cor ) {
          metXCor = metXInt;
          metYCor = metYInt;
          metReduc = (met2-met2Int)/met2Int; 
          met2Cor = met2Int;
          sumetCor = sumetInt;
          significanceCor = std::sqrt(met2Cor/sumetCor);
          iCor = i;
        }
      }
      //
      // If the MET must be significanly reduced, schedule the candidate to be cleaned
      if ( metReduc > minDeltaMet_ ) { 
        pfCandidatesToBeRemoved.push_back(iCor);
        metX = metXCor;
        metY = metYCor;
        sumet = sumetCor;
        met2 = met2Cor;
      } else { 
      // Otherwise just stop the loop
        next = false;
      }
    }
    //
    // The significance must be significantly reduced to indeed clean the candidates
    if ( significance - significanceCor > minSignificanceReduction_ && 
         significanceCor < maxSignificance_ ) {
      std::cout << "Significance reduction = " << significance << " -> " 
                << significanceCor << " = " << significanceCor - significance 
                << std::endl;
      for(unsigned j=0; j<pfCandidatesToBeRemoved.size(); ++j) {
        std::cout << "Removed : " << (*pfCandidates_)[pfCandidatesToBeRemoved[j]] << std::endl;
        pfCleanedCandidates_->push_back( (*pfCandidates_)[ pfCandidatesToBeRemoved[j] ] );
        (*pfCandidates_)[pfCandidatesToBeRemoved[j]].rescaleMomentum(1E-6);
        //reco::PFCandidate::ParticleType unknown = reco::PFCandidate::X;
        //(*pfCandidates_)[pfCandidatesToBeRemoved[j]].setParticleType(unknown);
      }
    }
  }

}

void 
PFAlgo::checkCleaning( const reco::PFRecHitCollection& cleanedHits ) { 
  

  // No hits to recover, leave.
  if ( !cleanedHits.size() ) return;

  //Compute met and met significance (met/sqrt(SumEt))
  double metX = 0.;
  double metY = 0.;
  double sumet = 0;
  std::vector<unsigned int> hitsToBeAdded;
  for(unsigned i=0; i<pfCandidates_->size(); i++) {
    const PFCandidate& pfc = (*pfCandidates_)[i];
    metX += pfc.px();
    metY += pfc.py();
    sumet += pfc.pt();
  }
  double met2 = metX*metX+metY*metY;
  double met2_Original = met2;
  // Select events with large MET significance.
  // double significance = std::sqrt(met2/sumet);
  // double significanceCor = significance;
  double metXCor = metX;
  double metYCor = metY;
  double sumetCor = sumet;
  double met2Cor = met2;
  bool next = true;
  unsigned iCor = 1E9;
  
  // Find the cleaned hit with the largest effect on the MET
  while ( next ) { 
    
    double metReduc = -1.;
    // Loop on the candidates
    for(unsigned i=0; i<cleanedHits.size(); ++i) {
      const PFRecHit& hit = cleanedHits[i];
      double length = std::sqrt(hit.position().Mag2()); 
      double px = hit.energy() * hit.position().x()/length;
      double py = hit.energy() * hit.position().y()/length;
      double pt = std::sqrt(px*px + py*py);
      
      // Check that it is  not already scheculed to be cleaned
      bool skip = false;
      for(unsigned j=0; j<hitsToBeAdded.size(); ++j) {
        if ( i == hitsToBeAdded[j] ) skip = true;
        if ( skip ) break;
      }
      if ( skip ) continue;
      
      // Now add the candidate to the MET
      double metXInt = metX + px;
      double metYInt = metY + py;
      double sumetInt = sumet + pt;
      double met2Int = metXInt*metXInt+metYInt*metYInt;

      // And check if it could contribute to a MET reduction
      if ( met2Int < met2Cor ) {
        metXCor = metXInt;
        metYCor = metYInt;
        metReduc = (met2-met2Int)/met2Int; 
        met2Cor = met2Int;
        sumetCor = sumetInt;
        // significanceCor = std::sqrt(met2Cor/sumetCor);
        iCor = i;
      }
    }
    //
    // If the MET must be significanly reduced, schedule the candidate to be added
    //
    if ( metReduc > minDeltaMet_ ) { 
      hitsToBeAdded.push_back(iCor);
      metX = metXCor;
      metY = metYCor;
      sumet = sumetCor;
      met2 = met2Cor;
    } else { 
      // Otherwise just stop the loop
      next = false;
    }
  }
  //
  // At least 10 GeV MET reduction
  if ( std::sqrt(met2_Original) - std::sqrt(met2) > 5. ) { 
    if ( debug_ ) { 
      std::cout << hitsToBeAdded.size() << " hits were re-added " << std::endl;
      std::cout << "MET reduction = " << std::sqrt(met2_Original) << " -> " 
                << std::sqrt(met2Cor) << " = " <<  std::sqrt(met2Cor) - std::sqrt(met2_Original)  
                << std::endl;
      std::cout << "Added after cleaning check : " << std::endl; 
    }
    for(unsigned j=0; j<hitsToBeAdded.size(); ++j) {
      const PFRecHit& hit = cleanedHits[hitsToBeAdded[j]];
      PFCluster cluster(hit.layer(), hit.energy(),
                        hit.position().x(), hit.position().y(), hit.position().z() );
      reconstructCluster(cluster,hit.energy());
      if ( debug_ ) { 
        std::cout << pfCandidates_->back() << ". time = " << hit.time() << std::endl;
      }
    }
  }

}


void 
PFAlgo::postMuonCleaning( const edm::Handle<reco::MuonCollection>& muonh,
                          const reco::VertexCollection& primaryVertices) {
  

  bool printout = false;
  // Check if the post muon Cleaning was requested - if not, do nothing
  //if ( !postMuonCleaning_ ) return;

  // Estimate SumEt from pile-up
  double sumetPU = 0.;
  for (unsigned short i=1 ;i<primaryVertices.size();++i ) {
    if ( !primaryVertices[i].isValid() || primaryVertices[i].isFake() ) continue; 
    primaryVertices[i];
    for ( reco::Vertex::trackRef_iterator itr = primaryVertices[i].tracks_begin();
          itr <  primaryVertices[i].tracks_end(); ++itr ) { 
      sumetPU += (*itr)->pt();
    }
  }
  sumetPU /= 0.65; // To estimate the neutral contribution from PU
  // std::cout << "Evaluation of sumetPU from " << primaryVertices.size() << " vertices : " << sumetPU << std::endl;

  //Compute met and met significance (met/sqrt(SumEt))
  double metX = 0.;
  double metY = 0.;
  double sumet = 0;
  double metXCosmics = 0.;
  double metYCosmics = 0.;
  double sumetCosmics = 0.;

  std::map<double,unsigned int> pfMuons;
  std::map<double,unsigned int> pfCosmics;
  typedef std::multimap<double, unsigned int>::iterator IE;

  for(unsigned i=0; i<pfCandidates_->size(); i++) {
    const PFCandidate& pfc = (*pfCandidates_)[i];
    double origin = 
      (pfc.vx()-primaryVertex_.x())*(pfc.vx()-primaryVertex_.x()) + 
      (pfc.vy()-primaryVertex_.y())*(pfc.vy()-primaryVertex_.y());

    // Compute MET
    metX += pfc.px();
    metY += pfc.py();
    sumet += pfc.pt();

    // Remember the muons and order them by decreasing energy
    // Remember the cosmic muons too 
    if ( pfc.particleId() == reco::PFCandidate::mu ) { 
      pfMuons.insert(std::pair<double,unsigned int>(-pfc.pt(),i));
      if ( origin > 1.0 ) { 
        pfCosmics.insert(std::pair<double,unsigned int>(-pfc.pt(),i));
        metXCosmics += pfc.px();
        metYCosmics += pfc.py();
        sumetCosmics += pfc.pt();
      }
    }
      
  }

  double met2 = metX*metX+metY*metY;
  double met2Cosmics = (metX-metXCosmics)*(metX-metXCosmics)+(metY-metYCosmics)*(metY-metYCosmics);

  // Clean cosmic muons
  if ( sumetCosmics > (sumet-sumetPU)/10. && met2Cosmics < met2 ) { 
    if ( printout ) 
      std::cout << "MEX,MEY,MET Before (Cosmics)" << metX << " " << metY << " " << std::sqrt(met2) << std::endl;
    for ( IE imu = pfCosmics.begin(); imu != pfCosmics.end(); ++imu ) { 
      const PFCandidate& pfc = (*pfCandidates_)[imu->second];
      //const PFCandidate pfcopy = pfc;
      pfCosmicsMuonCleanedCandidates_->push_back(pfc);
      if ( printout ) 
        std::cout << "Muon cleaned (cosmic). pt/d0 = " << pfc.pt() << " " 
                  << std::sqrt(pfc.vx()*pfc.vx() + pfc.vy()*pfc.vy()) << std::endl; 
      metX -= pfc.px();
      metY -= pfc.py();
      sumet -= pfc.pt();
      (*pfCandidates_)[imu->second].rescaleMomentum(1E-6);
    }
    met2 = metX*metX+metY*metY;
    if ( printout ) 
      std::cout << "MEX,MEY,MET After  (Cosmics)" << metX << " " << metY << " " << std::sqrt(met2) << std::endl;
  }

  // The remaining met
  // double metInitial = std::sqrt(met2);
  
  // Clean mismeasured muons
  for ( IE imu = pfMuons.begin(); imu != pfMuons.end(); ++imu ) { 
    const PFCandidate& pfc = (*pfCandidates_)[imu->second];

    // Don't care about low momentum muons
    if ( pfc.pt() < 20. ) continue;

    double metXNO = metX - pfc.px();
    double metYNO = metY - pfc.py();
    double met2NO = metXNO*metXNO + metYNO*metYNO; 
    double sumetNO = sumet - pfc.pt();
    double factor = std::max(2.,2000./sumetNO);

    reco::MuonRef muonRef = pfc.muonRef();
    reco::TrackRef bestMuTrack =  muonRef->muonBestTrack();
    if(!useBestMuonTrack_)
      bestMuTrack = muonRef->combinedMuon(); 
    reco::TrackRef trackerMu = muonRef->track();
    reco::TrackRef standAloneMu = muonRef->standAloneMuon();

    if ( !bestMuTrack || !trackerMu ) { 
      if ( sumetNO-sumetPU > 500. && met2NO < met2/4.) { 
        pfFakeMuonCleanedCandidates_->push_back(pfc);
        if ( printout ) {
          std::cout << "Muon cleaned (NO-bad) " << sumetNO << std::endl;
          std::cout << "sumet NO " << sumetNO << std::endl;
        }
        (*pfCandidates_)[imu->second].rescaleMomentum(1E-6);
        if ( printout ) 
          std::cout << "MEX,MEY,MET            " << metX << " " << metY << " " << std::sqrt(met2) << std::endl;
        metX = metXNO;
        metY = metYNO;
        met2 = met2NO;
        sumet = sumetNO;
        if ( printout ) {
          std::cout << "Muon cleaned (NO/TK) ! " << std::endl;
          std::cout << "MEX,MEY,MET Now (NO/TK)" << metX << " " << metY << " " << std::sqrt(met2) << std::endl;
        }
      }
    } else { 
      /*
      std::cout << "PF Muon vx,vy " << pfc.vx() << " " << pfc.vy() << std::endl;
      std::cout << "PF Muon px,py,pt: " << pfc.px() << " " << pfc.py() << " " << pfc.pt() << std::endl;
      std::cout << "TK Muon px,py,pt: " << trackerMu->px() << " " << trackerMu->py() << " " << trackerMu->pt() << " " << std::endl;
      std::cout << "GL Muon px,py,pt: " << bestMuTrack->px() << " " << bestMuTrack->py() << " " << bestMuTrack->pt() << " " << std::endl; 
      std::cout << "ST Muon px,py,pt: " << standAloneMu->px() << " " << standAloneMu->py() << " " << standAloneMu->pt() << " " << std::endl;
      std::cout << "isolated ? " << PFMuonAlgo::isIsolatedMuon(muonRef) << std::endl; 
      */
      double metXTK = metXNO + trackerMu->px(); 
      double metYTK = metYNO + trackerMu->py();
      double met2TK = metXTK*metXTK + metYTK*metYTK;
 
      double metXGL = metXNO + bestMuTrack->px();
      double metYGL = metYNO + bestMuTrack->py();
      double met2GL = metXGL*metXGL + metYGL*metYGL; 

      /*
      double metXST = metXNO + standAloneMu->px();
      double metYST = metYNO + standAloneMu->py();
      double met2ST = metXST*metXST + metYST*metYST; 
      */

      /*
      std::cout << "sumet NO " << sumetNO << std::endl;
      std::cout << "MEX,MEY,MET NO" << metXNO << " " << metYNO << " " << std::sqrt(met2NO) << std::endl;
      std::cout << "MEX,MEY,MET TK" << metXTK << " " << metYTK << " " << std::sqrt(met2TK) << std::endl;
      std::cout << "MEX,MEY,MET GL" << metXGL << " " << metYGL << " " << std::sqrt(met2GL) << std::endl;
      std::cout << "MEX,MEY,MET ST" << metXST << " " << metYST << " " << std::sqrt(met2ST) << std::endl;
      */

      //bool fixed = false;
      if ( ( sumetNO-sumetPU > 250. && met2TK < met2/4. && met2TK < met2GL ) || 
           ( met2TK < met2/2. && trackerMu->pt() < bestMuTrack->pt()/4. && met2TK < met2GL ) )  { 
        pfCleanedTrackerAndGlobalMuonCandidates_->push_back(pfc);
        math::XYZTLorentzVectorD p4(trackerMu->px(),trackerMu->py(),trackerMu->pz(),
                                    std::sqrt(trackerMu->p()*trackerMu->p()+0.1057*0.1057));
        if ( printout ) 
          std::cout << "Muon cleaned (TK) ! " << met2TK/met2 << " " << trackerMu->pt() / bestMuTrack->pt() << std::endl;
        (*pfCandidates_)[imu->second].setP4(p4);
        if ( printout ) 
          std::cout << "MEX,MEY,MET Before (TK)" << metX << " " << metY << " " << std::sqrt(met2) << std::endl;
        metX = metXTK;
        metY = metYTK;
        met2 = met2TK;
        //fixed = true;
        if ( printout ) 
          std::cout << "MEX,MEY,MET Now    (TK)" << metX << " " << metY << " " << std::sqrt(met2) << std::endl;
      } 
      
      else if ( ( sumetNO-sumetPU > 250. && met2GL < met2/4. && met2GL < met2TK ) ||  
                ( met2GL < met2/2. && bestMuTrack->pt() < trackerMu->pt()/4.&& met2GL < met2TK ) )  { 
        pfCleanedTrackerAndGlobalMuonCandidates_->push_back(pfc);
        math::XYZTLorentzVectorD p4(bestMuTrack->px(),bestMuTrack->py(),bestMuTrack->pz(),
                                    std::sqrt(bestMuTrack->p()*bestMuTrack->p()+0.1057*0.1057));
        if ( printout ) 
          std::cout << "Muon cleaned (GL) ! " << met2GL/met2 << " " << bestMuTrack->pt()/trackerMu->pt() <<  std::endl;
        (*pfCandidates_)[imu->second].setP4(p4);
        if ( printout ) 
          std::cout << "MEX,MEY,MET before (GL)" << metX << " " << metY << " " << std::sqrt(met2) << std::endl;
        metX = metXGL;
        metY = metYGL;
        met2 = met2GL;
        //fixed = true;
        if ( printout ) 
          std::cout << "MEX,MEY,MET Now    (GL)" << metX << " " << metY << " " << std::sqrt(met2) << std::endl;
      }

      // Fake muons at large pseudo-rapidity
      bool fake = 
        fabs ( pfc.eta() ) > 2.15 && 
        met2NO < met2/25. && 
        (met2GL < met2TK/2. || met2TK < met2GL/2.) &&
        standAloneMu->pt() < bestMuTrack->pt()/10. &&
        standAloneMu->pt() < trackerMu->pt()/10.; 

      // Fake and/or punch-through candidates
      bool punchthrough1 =  
        ( sumetNO-sumetPU > 250. && met2NO < met2/4. && (met2GL < met2TK/factor || met2TK < met2GL/factor) );

      // Now look for punch through candidates (i.e., muon with large neutral hadron behind)
      bool punchthrough2 =  
        pfc.p() > 100. &&  pfc.rawHcalEnergy() > 100. && 
        pfc.rawEcalEnergy()+pfc.rawHcalEnergy() > pfc.p()/3. &&
        !PFMuonAlgo::isIsolatedMuon(muonRef) && met2NO < met2/4.;

      if ( punchthrough1 ||  punchthrough2 || fake ) { 
        
        // Find the block of the muon
        const PFCandidate::ElementsInBlocks& eleInBlocks = pfc.elementsInBlocks();
        if ( !eleInBlocks.size() ) { 
          ostringstream err;
          err<<"A muon has no associated elements in block. Cannot proceed with postMuonCleaning()";
          edm::LogError("PFAlgo")<<err.str()<<endl;
        } else { 
          PFBlockRef blockRefMuon = eleInBlocks[0].first;
          unsigned indexMuon = eleInBlocks[0].second;
          for ( unsigned iele = 1; iele < eleInBlocks.size(); ++iele ) { 
            indexMuon = eleInBlocks[iele].second;
            break;
          }
          
          // Check if the muon gave rise to a neutral hadron
          double iHad = 1E9;
          bool hadron = false;
          for ( unsigned i = imu->second+1; i < pfCandidates_->size(); ++i ) { 
            const PFCandidate& pfcn = (*pfCandidates_)[i];
            const PFCandidate::ElementsInBlocks& ele = pfcn.elementsInBlocks();
            if ( !ele.size() ) { 
              ostringstream err2;
              err2<<"A pfCandidate has no associated elements in block. Cannot proceed with postMuonCleaning()";
              edm::LogError("PFAlgo")<<err2.str()<<endl;
              continue;
            }
            PFBlockRef blockRefHadron = ele[0].first;
            unsigned indexHadron = ele[0].second;
            // We are out of the block -> exit the loop
            if ( blockRefHadron.key() != blockRefMuon.key() ) break;
            // Check that this particle is a neutral hadron
            if ( indexHadron == indexMuon && 
                 pfcn.particleId() == reco::PFCandidate::h0 ) {
              iHad = i;
              hadron = true;
            }
            if ( hadron ) break;
          }
          
          if ( hadron ) { 
            const PFCandidate& pfch = (*pfCandidates_)[iHad];
            pfPunchThroughMuonCleanedCandidates_->push_back(pfc);
            pfPunchThroughHadronCleanedCandidates_->push_back(pfch);
            if ( printout ) {
              std::cout << pfc << std::endl;
              std::cout << "Raw ecal/hcal : " << pfc.rawEcalEnergy() << " " << pfc.rawHcalEnergy() << std::endl;
              std::cout << "Hadron: " << (*pfCandidates_)[iHad] << std::endl;
            }
            double rescaleFactor = (*pfCandidates_)[iHad].p()/(*pfCandidates_)[imu->second].p();
            metX -=  (*pfCandidates_)[imu->second].px() + (*pfCandidates_)[iHad].px();
            metY -=  (*pfCandidates_)[imu->second].py() + (*pfCandidates_)[iHad].py();
            (*pfCandidates_)[imu->second].rescaleMomentum(rescaleFactor);
            (*pfCandidates_)[iHad].rescaleMomentum(1E-6);
            (*pfCandidates_)[imu->second].setParticleType(reco::PFCandidate::h);
            if ( printout ) 
              std::cout << "MEX,MEY,MET Before " << metX << " " << metY << " " << std::sqrt(met2) << std::endl;
            metX +=  (*pfCandidates_)[imu->second].px();
            metY +=  (*pfCandidates_)[imu->second].py();          
            met2 = metX*metX + metY*metY;
            if ( printout ) {
              std::cout << "Muon changed to charged hadron" << std::endl;
              std::cout << "MEX,MEY,MET Now (NO)" << metX << " " << metY << " " << std::sqrt(met2) << std::endl;
            }
          } else if ( punchthrough1 || fake ) {
            const PFCandidate& pfc = (*pfCandidates_)[imu->second];
            pfFakeMuonCleanedCandidates_->push_back(pfc);
            if ( printout ) 
              std::cout << "MEX,MEY,MET Before " << metX << " " << metY << " " << std::sqrt(met2) << std::endl;
            metX -=  (*pfCandidates_)[imu->second].px();
            metY -=  (*pfCandidates_)[imu->second].py();          
            met2 = metX*metX + metY*metY;
            (*pfCandidates_)[imu->second].rescaleMomentum(1E-6);
            if ( printout ) {
              std::cout << "Muon cleaned (NO)" << std::endl;
              std::cout << "MEX,MEY,MET Now (NO)" << metX << " " << metY << " " << std::sqrt(met2) << std::endl;
            }
          }
        }
      }
      /*
      if ( !fixed ) { 
        std::cout << "TK Muon px,py,pt: " << trackerMu->px() << " " << trackerMu->py() << " " << trackerMu->pt() << " " << std::endl;
        std::cout << "GL Muon px,py,pt: " << bestMuTrack->px() << " " << bestMuTrack->py() << " " << bestMuTrack->pt() << " " << std::endl;
        std::cout << "MEX,MEY,MET PF " << metX << " " << metY << " " << std::sqrt(met2) << std::endl;
        std::cout << "MEX,MEY,MET NO " << metXNO << " " << metYNO << " " << std::sqrt(met2NO) << std::endl;
        std::cout << "MEX,MEY,MET TK " << metXTK << " " << metYTK << " " << std::sqrt(met2TK) << std::endl;
        std::cout << "MEX,MEY,MET GL " << metXGL << " " << metYGL << " " << std::sqrt(met2GL) << std::endl;
      }
      */
    }

  }

  // And now, add muons which did not make it for various reasons
  for ( unsigned imu = 0; imu < muonh->size(); ++imu ) {
    reco::MuonRef muonRef( muonh, imu );
    // If not used, check its effect on met
    reco::TrackRef bestMuTrack =  muonRef->muonBestTrack();
    if(!useBestMuonTrack_)
      bestMuTrack = muonRef->combinedMuon(); 
    reco::TrackRef trackerMu = muonRef->track();
    reco::TrackRef standAloneMu = muonRef->standAloneMuon();

    // check if the muons has already been taken
    bool used = false;
    bool hadron = false;
    for(unsigned i=0; i<pfCandidates_->size(); i++) {
      const PFCandidate& pfc = (*pfCandidates_)[i];
      if ( !pfc.trackRef().isNonnull() ) continue;

      if ( pfc.trackRef().isNonnull() && pfc.trackRef() == trackerMu ) { 
        hadron = true;
      }

      // The pf candidate is not associated to a muon
      // if ( !pfc.muonRef().isNonnull() || pfc.muonRef() != muonRef ) continue; ! This test is buggy !     
      if ( !pfc.muonRef().isNonnull() ) continue;
      
      // Check if the muon is used... 
      if ( pfc.muonRef()->track() == trackerMu || pfc.muonRef()->combinedMuon() == bestMuTrack ) {
        if ( printout ) { 
          std::cout << "This muon is already used ..." << std::endl;
          std::cout << pfc << std::endl;
          std::cout << muonRef->p() << " " << muonRef->pt() << " " << muonRef->eta() << " " << muonRef->phi() << std::endl;
        }
        used = true;
      } 
      else {
        // Check if the stand-alone muon is not a spurious copy of an existing muon 
        // (Protection needed for HLT)
        if ( pfc.muonRef()->isStandAloneMuon() && muonRef->isStandAloneMuon() ) { 
          double dEta = pfc.muonRef()->standAloneMuon()->eta() - standAloneMu->eta();
          double dPhi = pfc.muonRef()->standAloneMuon()->phi() - standAloneMu->phi();
          double dR = sqrt(dEta*dEta + dPhi*dPhi);
          if ( printout ) { 
            std::cout << "StandAlone to be added ? " << std::endl;
            std::cout << " eta = " << pfc.muonRef()->standAloneMuon()->eta() << " " << standAloneMu->eta() << std::endl;
            std::cout << " phi = " << pfc.muonRef()->standAloneMuon()->phi() << " " << standAloneMu->phi() << std::endl;
            std::cout << " pt = " << pfc.muonRef()->standAloneMuon()->pt() << " " << standAloneMu->pt() << std::endl;
            std::cout << " Delta R = " << dR << std::endl;
          }
          if ( dR < 0.005 ) { 
            used = true;
            if ( printout ) std::cout << "Not removed !" << std::endl;
          }
        }
      } 
      if ( used ) break; 

    }

    if ( used ) continue;

    double ptGL = muonRef->isGlobalMuon() ? bestMuTrack->pt() : 0.;
    double pxGL = muonRef->isGlobalMuon() ? bestMuTrack->px() : 0.;
    double pyGL = muonRef->isGlobalMuon() ? bestMuTrack->py() : 0.;
    double pzGL = muonRef->isGlobalMuon() ? bestMuTrack->pz() : 0.;

    double ptTK = muonRef->isTrackerMuon() ? trackerMu->pt() : 0.;
    double pxTK = muonRef->isTrackerMuon() ? trackerMu->px() : 0.;
    double pyTK = muonRef->isTrackerMuon() ? trackerMu->py() : 0.;
    double pzTK = muonRef->isTrackerMuon() ? trackerMu->pz() : 0.;

    double ptST = muonRef->isStandAloneMuon() ? standAloneMu->pt() : 0.;
    double pxST = muonRef->isStandAloneMuon() ? standAloneMu->px() : 0.;
    double pyST = muonRef->isStandAloneMuon() ? standAloneMu->py() : 0.;
    double pzST = muonRef->isStandAloneMuon() ? standAloneMu->pz() : 0.;

    //std::cout << "pT TK/GL/ST : " << ptTK << " " << ptGL << " " << ptST << std::endl;

    double metXTK = metX + pxTK; 
    double metYTK = metY + pyTK;
    double met2TK = metXTK*metXTK + metYTK*metYTK;
    
    double metXGL = metX + pxGL;
    double metYGL = metY + pyGL;
    double met2GL = metXGL*metXGL + metYGL*metYGL; 

    double metXST = metX + pxST; 
    double metYST = metY + pyST;
    double met2ST = metXST*metXST + metYST*metYST;

    //std::cout << "met TK/GL/ST : " << sqrt(met2) << " " << sqrt(met2TK) << " " << sqrt(met2GL) << " " << sqrt(met2ST) << std::endl;


    if ( ptTK > 20. && met2TK < met2/4. && met2TK < met2GL && met2TK < met2ST ) { 
      double energy = std::sqrt(pxTK*pxTK+pyTK*pyTK+pzTK*pzTK+0.1057*0.1057);
      int charge = trackerMu->charge()>0 ? 1 : -1;
      math::XYZTLorentzVector momentum(pxTK,pyTK,pzTK,energy);
      reco::PFCandidate::ParticleType particleType = 
         ptGL > 20. ? reco::PFCandidate::mu : reco::PFCandidate::h;
      double radius = std::sqrt(
         (trackerMu->vertex().x()-primaryVertex_.x())*(trackerMu->vertex().x()-primaryVertex_.x())+
         (trackerMu->vertex().y()-primaryVertex_.y())*(trackerMu->vertex().y()-primaryVertex_.y()));
      
      // Add it to the stack ( if not already in the list ) 
      if ( !hadron && radius < 1.0 ) { 
        pfCandidates_->push_back( PFCandidate( charge, 
                                               momentum,
                                               particleType ) );
        
        pfCandidates_->back().setVertexSource( PFCandidate::kTrkMuonVertex );
        pfCandidates_->back().setTrackRef( trackerMu );
        pfCandidates_->back().setMuonRef( muonRef );
        
        if ( printout ) {
          std::cout << "MEX,MEY,MET before " << metX << " " << metY << " " << std::sqrt(met2) << std::endl;
          std::cout << "Muon TK added " << std::endl;
          std::cout << "pT TK/GL/ST : " << ptTK << " " << ptGL << " " << ptST << std::endl;
        }
        metX +=  pfCandidates_->back().px();
        metY +=  pfCandidates_->back().py();      
        met2 = metX*metX + metY*metY;
        if ( printout ) {
          std::cout <<  pfCandidates_->back() << std::endl;
          std::cout << "MEX,MEY,MET   now " << metX << " " << metY << " " << std::sqrt(met2) << std::endl; 
        }

        const PFCandidate& pfc = pfCandidates_->back();
        pfAddedMuonCandidates_->push_back(pfc);

      }

    } else  if (  ptGL > 20. && met2GL < met2/4. && met2GL < met2TK && met2GL < met2ST ) { 

      double energy = std::sqrt(pxGL*pxGL+pyGL*pyGL+pzGL*pzGL+0.1057*0.1057);
      int charge = bestMuTrack->charge()>0 ? 1 : -1;
      math::XYZTLorentzVector momentum(pxGL,pyGL,pzGL,energy);
      reco::PFCandidate::ParticleType particleType = reco::PFCandidate::mu;
      double radius = std::sqrt(
        (bestMuTrack->vertex().x()-primaryVertex_.x())*(bestMuTrack->vertex().x()-primaryVertex_.x())+
        (bestMuTrack->vertex().y()-primaryVertex_.y())*(bestMuTrack->vertex().y()-primaryVertex_.y()));
      
      // Add it to the stack
      if ( radius < 1.0 ) { 
        pfCandidates_->push_back( PFCandidate( charge, 
                                               momentum,
                                               particleType ) );
        
        pfCandidates_->back().setVertexSource( PFCandidate::kComMuonVertex );
        //if ( ptTK > 0. ) 
        if (trackerMu.isNonnull() ) pfCandidates_->back().setTrackRef( trackerMu );
        pfCandidates_->back().setMuonRef( muonRef );
        
        if ( printout ) {
          std::cout << "MEX,MEY,MET before " << metX << " " << metY << " " << std::sqrt(met2) << std::endl;
          std::cout << "Muon GL added " << std::endl;
          std::cout << "pT TK/GL/ST : " << ptTK << " " << ptGL << " " << ptST << std::endl;
        }
        
        metX +=  pfCandidates_->back().px();
        metY +=  pfCandidates_->back().py();      
        met2 = metX*metX + metY*metY;
        
        if ( printout ) {
          std::cout <<  pfCandidates_->back() << std::endl;
          std::cout << "MEX,MEY,MET   now " << metX << " " << metY << " " << std::sqrt(met2) << std::endl;
        }

        const PFCandidate& pfc = pfCandidates_->back();
        pfAddedMuonCandidates_->push_back(pfc);

      }
        
    } else  if ( ptST > 20. &&  met2ST < met2/4. && met2ST < met2TK && met2ST < met2GL ) { 

      double energy = std::sqrt(pxST*pxST+pyST*pyST+pzST*pzST+0.1057*0.1057);
      int charge = standAloneMu->charge()>0 ? 1 : -1;
      math::XYZTLorentzVector momentum(pxST,pyST,pzST,energy);
      reco::PFCandidate::ParticleType particleType = reco::PFCandidate::mu;
      double radius = std::sqrt(
        (standAloneMu->vertex().x()-primaryVertex_.x())*(standAloneMu->vertex().x()-primaryVertex_.x())+
        (standAloneMu->vertex().y()-primaryVertex_.y())*(standAloneMu->vertex().y()-primaryVertex_.y()));
      
      // Add it to the stack
      if ( radius < 1.0 ) { 
        pfCandidates_->push_back( PFCandidate( charge, 
                                               momentum,
                                               particleType ) );
        
        pfCandidates_->back().setVertexSource( PFCandidate::kSAMuonVertex);
        if (trackerMu.isNonnull() ) pfCandidates_->back().setTrackRef( trackerMu );
        pfCandidates_->back().setMuonRef( muonRef );
        
        if ( printout ) {
          std::cout << "MEX,MEY,MET before " << metX << " " << metY << " " << std::sqrt(met2) << std::endl;
          std::cout << "Muon ST added " << std::endl;
          std::cout << "pT TK/GL/ST : " << ptTK << " " << ptGL << " " << ptST << std::endl;
        }

        metX +=  pfCandidates_->back().px();
        metY +=  pfCandidates_->back().py();      
        met2 = metX*metX + metY*metY;
        
        if ( printout ) {
          std::cout <<  pfCandidates_->back() << std::endl;
          std::cout << "MEX,MEY,MET   now " << metX << " " << metY << " " << std::sqrt(met2) << std::endl;
        }

        const PFCandidate& pfc = pfCandidates_->back();
        pfAddedMuonCandidates_->push_back(pfc);

      }
    }

  }

  /*
  if ( std::sqrt(met2) > 500. ) { 
    std::cout << "MET initial : " << metInitial << std::endl;
    std::cout << "MET final : " << std::sqrt(met2) << std::endl;
  }
  */

}

void PFAlgo::setElectronExtraRef(const edm::OrphanHandle<reco::PFCandidateElectronExtraCollection >& extrah) {
  if(!usePFElectrons_) return;
  //  std::cout << " setElectronExtraRef " << std::endl;
  unsigned size=pfCandidates_->size();

  for(unsigned ic=0;ic<size;++ic) {
    // select the electrons and add the extra
    if((*pfCandidates_)[ic].particleId()==PFCandidate::e) {
      
      PFElectronExtraEqual myExtraEqual((*pfCandidates_)[ic].gsfTrackRef());
      std::vector<PFCandidateElectronExtra>::const_iterator it=find_if(pfElectronExtra_.begin(),pfElectronExtra_.end(),myExtraEqual);
      if(it!=pfElectronExtra_.end()) {
        //      std::cout << " Index " << it-pfElectronExtra_.begin() << std::endl;
        reco::PFCandidateElectronExtraRef theRef(extrah,it-pfElectronExtra_.begin());
        (*pfCandidates_)[ic].setPFElectronExtraRef(theRef);
      }
      else {
        (*pfCandidates_)[ic].setPFElectronExtraRef(PFCandidateElectronExtraRef());
      }
    }
    else  // else save the mva and the extra as well ! 
      {
        if((*pfCandidates_)[ic].trackRef().isNonnull()) {
          PFElectronExtraKfEqual myExtraEqual((*pfCandidates_)[ic].trackRef());
          std::vector<PFCandidateElectronExtra>::const_iterator it=find_if(pfElectronExtra_.begin(),pfElectronExtra_.end(),myExtraEqual);
          if(it!=pfElectronExtra_.end()) {
            (*pfCandidates_)[ic].set_mva_e_pi(it->mvaVariable(PFCandidateElectronExtra::MVA_MVA));
            reco::PFCandidateElectronExtraRef theRef(extrah,it-pfElectronExtra_.begin());
            (*pfCandidates_)[ic].setPFElectronExtraRef(theRef);
            (*pfCandidates_)[ic].setGsfTrackRef(it->gsfTrackRef());
          }     
        }
      }

  }

  size=pfElectronCandidates_->size();
  for(unsigned ic=0;ic<size;++ic) {
    // select the electrons - this test is actually not needed for this collection
    if((*pfElectronCandidates_)[ic].particleId()==PFCandidate::e) {
      // find the corresponding extra
      PFElectronExtraEqual myExtraEqual((*pfElectronCandidates_)[ic].gsfTrackRef());
      std::vector<PFCandidateElectronExtra>::const_iterator it=find_if(pfElectronExtra_.begin(),pfElectronExtra_.end(),myExtraEqual);
      if(it!=pfElectronExtra_.end()) {
        reco::PFCandidateElectronExtraRef theRef(extrah,it-pfElectronExtra_.begin());
        (*pfElectronCandidates_)[ic].setPFElectronExtraRef(theRef);

      }
    }
  }

}
void PFAlgo::setPhotonExtraRef(const edm::OrphanHandle<reco::PFCandidatePhotonExtraCollection >& ph_extrah) {
  if(!usePFPhotons_) return;  
  unsigned int size=pfCandidates_->size();
  unsigned int sizePhExtra = pfPhotonExtra_.size();
  for(unsigned ic=0;ic<size;++ic) {
    // select the electrons and add the extra
    if((*pfCandidates_)[ic].particleId()==PFCandidate::gamma && (*pfCandidates_)[ic].mva_nothing_gamma() > 0.99) {
      if((*pfCandidates_)[ic].superClusterRef().isNonnull()) {
        bool found = false;
        for(unsigned pcextra=0;pcextra<sizePhExtra;++pcextra) {
          if((*pfCandidates_)[ic].superClusterRef() == pfPhotonExtra_[pcextra].superClusterRef()) {
            reco::PFCandidatePhotonExtraRef theRef(ph_extrah,pcextra);
            (*pfCandidates_)[ic].setPFPhotonExtraRef(theRef);
            found = true;
            break;
          }
        }
        if(!found) 
          (*pfCandidates_)[ic].setPFPhotonExtraRef((PFCandidatePhotonExtraRef())); // null ref
      }
      else {
        (*pfCandidates_)[ic].setPFPhotonExtraRef((PFCandidatePhotonExtraRef())); // null ref
      }
    }
  }      
}