---

PFAlgo.cc

Go to the documentation of this file.

#include "RecoParticleFlow/PFAlgo/interface/PFAlgo.h"
#include "RecoParticleFlow/PFAlgo/interface/PFElectronAlgo.h"  //PFElectrons
#include "RecoParticleFlow/PFAlgo/interface/PFConversionAlgo.h"  // PFConversions

#include "RecoParticleFlow/PFClusterTools/interface/PFEnergyCalibration.h"
#include "RecoParticleFlow/PFClusterTools/interface/PFClusterCalibration.h" // Brand new cluster calibration !

#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/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/ParticleFlowCandidate/interface/PFCandidate.h"
#include "DataFormats/ParticleFlowCandidate/interface/PFCandidateFwd.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;

float PFAlgo::maxMvaCut_ = PFCandidate::bigMva_;


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



PFAlgo::PFAlgo()
  : pfCandidates_( new PFCandidateCollection),
    nSigmaECAL_(3), 
    nSigmaHCAL_(3),
    clusterRecovery_(false), 
    PSCut_(999.),
    mergedPhotonsMVA_( 0 ),
    mvaCut_(PFAlgo::maxMvaCut_),
    algo_(1),
    debug_(false),
    pfele_(0)
{}

PFAlgo::~PFAlgo() {
  if (usePFElectrons_) delete pfele_;
  if (usePFConversions_) delete pfConversion_;
}


void 
PFAlgo::setParameters(double nSigmaECAL,
                      double nSigmaHCAL, 
                      const shared_ptr<PFEnergyCalibration>& calibration,
                      const shared_ptr<pftools::PFClusterCalibration>& clusterCalibration,
                      bool   newCalib,
                      bool   clusterRecovery,
                      double PSCut,
                      double mvaCut,
                      const  char* mvaWeightFile ) {

  nSigmaECAL_ = nSigmaECAL;
  nSigmaHCAL_ = nSigmaHCAL;

  calibration_ = calibration;
  clusterCalibration_ = clusterCalibration;
  newCalib_ = newCalib;
  
  // std::cout << "Cluster calibration parameters : " << *clusterCalibration_ << std::endl;

  clusterRecovery_ = clusterRecovery;

  PSCut_=PSCut;
  mvaCut_ = mvaCut;
  
  if( mergedPhotonsMVA_ ) delete mergedPhotonsMVA_;
  mergedPhotonsMVA_ = new TMVA::Reader();
  mergedPhotonsMVA_->AddVariable("eECAL/pTrack", &eECALOverpTrack_);
  mergedPhotonsMVA_->AddVariable("chi2ECAL", &chi2ECAL_);
  mergedPhotonsMVA_->AddVariable("ptTrack", &ptTrack_);
  
  if( mvaCut < PFAlgo::maxMvaCut_ ) {
    FILE * file = fopen(mvaWeightFile, "r");
    if (file) {
      fclose(file);
      // file is readable, book MVA
      mergedPhotonsMVA_->BookMVA( "MVA", mvaWeightFile );
    }
    else {
      // weight file is not readable
      string err = "PFAlgo: cannot open weight file '";
      err += mvaWeightFile;
      err += "'";
      throw invalid_argument( err );
    }
  }
}

//PFElectrons: a new method added to set the parameters for electron reconstruction. 
void 
PFAlgo::setPFEleParameters(double chi2EcalGSF,
                           double chi2EcalBrem,
                           double chi2HcalGSF,
                           double chi2HcalBrem,
                           double chi2PsGSF,
                           double chi2PsBrem,
                           double mvaEleCut,
                           string mvaWeightFileEleID,
                           bool usePFElectrons) {
  setchi2Values_.push_back(chi2EcalGSF);
  setchi2Values_.push_back(chi2EcalBrem);
  setchi2Values_.push_back(chi2HcalGSF);
  setchi2Values_.push_back(chi2HcalBrem);
  setchi2Values_.push_back(chi2PsGSF);
  setchi2Values_.push_back(chi2PsBrem);

  mvaEleCut_ = mvaEleCut;
  usePFElectrons_ = usePFElectrons;
  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(setchi2Values_,mvaEleCut_,mvaWeightFileEleID_);
}


void PFAlgo::setPFConversionParameters(bool usePFConversions ) {

  usePFConversions_ = usePFConversions;
  pfConversion_ = new PFConversionAlgo();


}



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( 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 > 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(!singleEcalOrHcal) {
      otherBlockRefs.push_back( blockref );
    }
  }//loop blocks
  
  if( debug_ ){
    cout<<"# Ecal blocks: "<<ecalBlockRefs.size()
        <<", # Hcal blocks: "<<hcalBlockRefs.size()
        <<", # Other blocks: "<<otherBlockRefs.size()<<endl;
  }


  // loop on blocks that are not single ecal, 
  // and not single hcal.
  // single ecal and single hcal will be used in cluster recovery

  for( IBR io = otherBlockRefs.begin(); io!=otherBlockRefs.end(); ++io) {
    processBlock( *io, hcalBlockRefs, ecalBlockRefs );
  }

  // cluster recovery has been done. 
  // process remaining ecal and hcal clusters. 
  
  std::list< reco::PFBlockRef > empty;

  // process remaining single hcal blocks
  for( IBR ih = hcalBlockRefs.begin(); ih!=hcalBlockRefs.end(); ++ih) {
    processBlock( *ih, empty, empty );
  }

  // process remaining single ecal blocks
  for( IBR ie = ecalBlockRefs.begin(); ie!=ecalBlockRefs.end(); ++ie) {
    processBlock( *ie, empty, empty );
  }
}


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

  typedef std::multimap<double, unsigned>::iterator IE;

  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
  if (usePFElectrons_) {
    if (pfele_->isElectronValidCandidate(blockref,active)){
      // if there is at least a valid candidate is get the vector of pfcandidates
      std::vector<reco::PFCandidate> PFElectCandidates_;
      PFElectCandidates_ = pfele_->getElectronCandidates();
      for ( std::vector<reco::PFCandidate>::iterator ec = PFElectCandidates_.begin();
            ec != PFElectCandidates_.end(); ++ec )
        {
          pfCandidates_->push_back(*ec);
          // the pfalgo candidates vector is filled
        }
      // The vector active is automatically changed (it is passed by ref) in PFElectronAlgo
      // for all the electron candidate
    }
  }

  
  // usePFConversions_ is used to switch ON/OFF the use of the PFConversionAlgo
  if (usePFConversions_) {
    if (pfConversion_->isConversionValidCandidate(blockref, active )){
      // if at least one conversion candidate is found ,it is fed to the final list of Pflow Candidates 
      std::vector<reco::PFCandidate> PFConversionCandidates = pfConversion_->conversionCandidates();
            
      for ( std::vector<reco::PFCandidate>::iterator iConv = PFConversionCandidates.begin(); 
            iConv != PFConversionCandidates.end(); ++iConv ) {
        pfCandidates_->push_back(*iConv);
      }
      for(unsigned iEle=0; iEle<elements.size(); iEle++) {
        //cout << " PFAlgo::processBlock conversion element " << iEle << " active " << active[iEle] << endl;
      }
    }
  }



  if(debug_) 
    cout<<endl<<"--------------- loop 1 ---------------------"<<endl;

  // 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 hcal and ecal elements 
  vector<unsigned> hcalIs;
  vector<unsigned> ecalIs;
  vector<unsigned> trackIs;


  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;
    default:
      continue;
    }

    // we're now dealing with a track
    unsigned iTrack = iEle;
    if(debug_) { 
      cout<<"TRACK"<<endl;
      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() ); 
                      

    // look for associated elements of all types

    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 ); 
    //MICHELE 
    //TEMPORARY SOLUTION FOR ELECTRON REJECTION IN PFTAU
    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;
    }

    typedef std::multimap<double, unsigned>::iterator IE;
   
    // loop on associated elements

    bool ecalFound = false;
    bool hcalFound = false;

    if(debug_) 
      cout<<"now looping on elements associated to the track"<<endl;
          
    // for each element associated to track iTrack

    double ecalEnergy = 0;

    unsigned iEcal = 99999;
    bool     connectedToEcal = false;

    for(IE ie = ecalElems.begin(); ie != ecalElems.end(); ++ie ) {
      
      // double  dist  = ie->first;
      unsigned index = ie->second;
      
      PFBlockElement::Type type = elements[index].type();
      
      if(debug_) {
        double chi2 = block.chi2(iTrack,index,linkData,reco::PFBlock::LINKTEST_ALL);
        cout<<"\telement "<<elements[index]<<" linked with chi2 "<<chi2<<endl;
      }

      assert( type == PFBlockElement::ECAL );

      // PFElectrons:Added a check on the activate elements
      // if the ECAL cluster has already been locked break the loop?y
      if ( ! active[index] ) break;      
      
      reco::PFClusterRef clusterRef = elements[index].clusterRef();
      assert( !clusterRef.isNull() );

      if( !ecalFound ) { 
        
        if(debug_) 
          cout<<"\t\tfirst ecal element "<<endl;
        
        ecalFound = true;
        
        if( hcalElems.empty() ) { 
          
          if(debug_) {
            cout<<"\t\tlocking ecal cluster"<<endl;
          }
            
          ecalEnergy    = clusterRef->energy();
          iEcal = index;
          connectedToEcal = true;

          // colin  in fact the cluster should maybe not be locked,
          // but some of its energy could make it to the hcal 
          // loop
          // PJ Colin is right - this is dealt with a little later.
          // active[index] = false;

        } //  if( hcalElems.empty() { 
        else {
          if(debug_) 
            cout<<"\t\tat least one hcal element connected to the track."
                <<" Sparing Ecal cluster for the hcal loop"<<endl; 
        }
      }
      else { // other associated ecal
        // this should not happen, except for very large 
        // chi2 max for track-ecal association
        // indeed, ecal clusters can't be too close
        // since seed look-up is done on 8 neighbours
        
        // anyway, when this happens, we cannot 
        // form a photon from the cluster
        // since this cluster might be connected to another 
        // track, or to an hcal cluster. 
        // so we do nothing.
        
        if(debug_) 
          cout<<"\t\tsecondary ecal element "<<endl;
      }
    
    }//loop ecal elements

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

    if( hcalElems.empty() ) {

      vector<unsigned> elementIndices;
      elementIndices.push_back(iTrack); 

      //  create a track candidate 
      
      active[iTrack] = false;
      
      // PJ If there is no HCAL, just remove the link by rechit with this track
      // Just trust the ECAL granularity
    
      
      unsigned tmpi = reconstructTrack( elements[iTrack] );
      
      double calibEcal = ecalEnergy;
      double calibHcal = 0.;
      double slopeEcal = 1.;
      if ( connectedToEcal && active[iEcal] ) {
        if ( newCalib_ ) { 
          clusterCalibration_->
            getCalibratedEnergyEmbedAInHcal(calibEcal, calibHcal,
                                            elements[iEcal].clusterRef()->positionREP().Eta(),
                                            elements[iEcal].clusterRef()->positionREP().Phi());
          if ( ecalEnergy > 0. ) slopeEcal = calibEcal/ecalEnergy;
        } else { 
          slopeEcal = calibration_->paramECALplusHCAL_slopeECAL();
          calibEcal *= slopeEcal; 
        }
      }

      (*pfCandidates_)[tmpi].setEcalEnergy( calibEcal );
      (*pfCandidates_)[tmpi].setHcalEnergy( 0 );
      (*pfCandidates_)[tmpi].setPs1Energy( 0 );
      (*pfCandidates_)[tmpi].setPs2Energy( 0 );
      (*pfCandidates_)[tmpi].addElementInBlock( blockref, iTrack ); 

      // Create some photon/neutral hadron (?) if the ecal energy is too large
      if( connectedToEcal && active[iEcal] ) {
        reco::PFClusterRef pivotalRef = elements[iEcal].clusterRef();
        (*pfCandidates_)[tmpi].addElementInBlock( blockref, iEcal ); 

        std::multimap<double, unsigned> sortedTracks;
        block.associatedElements( iEcal,  linkData,
                                  sortedTracks,
                                  reco::PFBlockElement::TRACK,
                                  reco::PFBlock::LINKTEST_ALL );
        // Check if the track closest to this ECAL is the current track
        // otherwise keeps this ECAL for later
        if ( sortedTracks.begin()->second == iTrack ) {  
          double neutralEnergy = calibEcal;

          // Check if there are other tracks linked to that ECAL
          for(IE ie = sortedTracks.begin(); ie != sortedTracks.end(); ++ie ) {
            unsigned jTrack = ie->second;
            
            // 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
          
          // PJ Don't forget to unfold the ECAL calibration : 01/10/08 !!!
          active[iEcal] = false;
          neutralEnergy /= slopeEcal;
          if ( neutralEnergy > 0. ) { 
            unsigned tmpj = reconstructCluster( *pivotalRef, neutralEnergy ); 
            (*pfCandidates_)[tmpj].setEcalEnergy( neutralEnergy );
            (*pfCandidates_)[tmpj].setHcalEnergy( 0. );
            (*pfCandidates_)[tmpj].setPs1Energy( -1 );
            (*pfCandidates_)[tmpj].setPs2Energy( -1 );
            (*pfCandidates_)[tmpj].addElementInBlock(blockref, iEcal);
            (*pfCandidates_)[tmpj].addElementInBlock( blockref, iTrack ); 
          } // End neutral energy
        } // End closest track
      } // End connected ECAL
      
      //MICHELE 
      //TEMPORARY SOLUTION FOR ELECTRON REJECTION IN PFTAU
       if (usePFElectrons_ == false) (*pfCandidates_)[tmpi].set_mva_e_pi(gsfElems.size()>0);
      //

    }

    for(IE ie = hcalElems.begin(); ie != hcalElems.end(); ++ie ) {
      
      // double   chi2  = ie->first;
      unsigned index = ie->second;
      
      PFBlockElement::Type type = elements[index].type();

      if(debug_) {
        double chi2 = block.chi2(iTrack,index,linkData,reco::PFBlock::LINKTEST_ALL);
        cout<<"\telement "<<elements[index]<<" linked with chi2 "<<chi2<<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, -1, linkData );
        //Note ALex: a fonction should be implemented to
        //force link to -1 for all possible tests
        //For the moment (15/11/2007) only 2 possible
        //tests: CHI2 and RECHIT.
        block.setLink( iTrack, index,  -1, -1, linkData, 
                       PFBlock::LINKTEST_CHI2 );
        block.setLink( iTrack, index,  -1, -1, linkData,
                       PFBlock::LINKTEST_RECHIT );
      }
    } //loop hcal elements   
  } // end of loop 1 on elements iEle of any type



  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;

    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
    // This printout pops up when a cluster was recovered by a previous track
    // It shows that the cluster recovery must be entirely revisited, 
    // and should be de-activated for now.
    //
    // if ( !active[iHcal] ) std::cout << "Warning !!! HCAL locked !!!" << std::endl;
    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 totalEcal = -0.0000000001; // thanks alex
    double totalHcal = hclusterref->energy();
    vector<double> hcalP;
    vector<double> hcalDP;
    vector<unsigned> tkIs;
    double maxDPovP = -9999.;
    
    //saving reference to associated tracks
    std::vector< reco::PFRecTrackRef > associatedTracks;
    
    vector< unsigned > chargedHadronsIndices;
    double mergedNeutralHadronEnergy = 0;

    for(IE ie = sortedTracks.begin(); ie != sortedTracks.end(); ++ie ) {
      // double chi2 = ie->first;


      unsigned iTrack = ie->second;
      //       tkIs.push_back(iTrack);

      // Check the track has not already been used (e.g., in electrons, conversions...)
      if ( !active[iTrack] ) continue;
           
      PFBlockElement::Type type = elements[iTrack].type();
      assert( type == reco::PFBlockElement::TRACK );
      
      reco::TrackRef trackRef = elements[iTrack].trackRef();
      assert( !trackRef.isNull() ); 
      
      unsigned tmpi = reconstructTrack( elements[iTrack] );
      
      reco::PFCandidate& currentChargedHadron = (*pfCandidates_)[tmpi];
//       currentChargedHadron.addElement(&elements[iTrack]);
//       currentChargedHadron.addElement(&elements[iHcal]);
      currentChargedHadron.addElementInBlock( blockref, iTrack );
      currentChargedHadron.addElementInBlock( blockref, iHcal );

      chargedHadronsIndices.push_back( tmpi );
      //MICHELE 
      //TEMPORARY SOLUTION FOR ELECTRON REJECTION IN PFTAU
      if (usePFElectrons_ == false) {
        std::multimap<double,unsigned> gsfElems;
        block.associatedElements( iTrack,  linkData,
                                  gsfElems ,
                                  reco::PFBlockElement::GSF );
        currentChargedHadron.set_mva_e_pi(gsfElems.size()>0);
      }
      //

      //filling reference to track vector
      //for cluster recovery algorithm
      associatedTracks.push_back( elements[iTrack].trackRefPF() );

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

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

      double Dp = trackRef->qoverpError()*trackMomentum*trackMomentum;
      sumpError2 += Dp*Dp;
      hcalP.push_back(trackMomentum);
      hcalDP.push_back(Dp);
      if (Dp/trackMomentum > maxDPovP) maxDPovP = Dp/trackMomentum;

      // look for ecal elements associated to iTrack (associated to iHcal)
      // PJ Possibility of link by rechit/chi2
      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;
      
      if( !sortedEcals.empty() ) 
        { // start case: at least one ecal element associated to iTrack
          // double chi2Ecal = sortedEcals.begin()->first;
          
          // PJ Check if this ECAL is not closer to another track
          // PJ If so take the next one, etc.
          for ( IE iec=sortedEcals.begin(); iec!=sortedEcals.end(); ++iec ) {

            unsigned iEcal = iec->second;
            
            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);
            if(debug_) cout<<"\t\t\tclosest: "
                           <<elements[iEcal]<<endl;
            
            if( active[iEcal] ) { // start case: iEcal is active
              PFBlockElement::Type type = elements[ iEcal ].type();
              assert( type == PFBlockElement::ECAL );
              PFClusterRef eclusterref = elements[iEcal].clusterRef();
              assert(!eclusterref.isNull() );
              
              
              totalEcal += eclusterref->energy();
              
              // ecal part of abc calibrated energy.
              // the offset is fully given to hcal energy
              float ecalEnergyCalibrated = eclusterref->energy();
              ecalEnergyCalibrated*=calibration_->paramECALplusHCAL_slopeECAL(); 
              
              currentChargedHadron.setEcalEnergy( ecalEnergyCalibrated );
              //             currentChargedHadron.addElement(&elements[iEcal]);
              currentChargedHadron.addElementInBlock( blockref, iEcal );
              
              std::pair<double, unsigned> associatedEcal 
                = make_pair( chi2Ecal, iEcal );
              associatedEcals.insert( make_pair(iTrack, associatedEcal) );
              active[iEcal] = false;
              
              if(debug_) cout<<"\t\t\tactive, adding "<<eclusterref->energy()
                             <<" to ECAL energy, and locking"<<endl;

            } // end case : iEcal is active
            else {
              if(debug_) cout<<"\t\tcluster locked"<<endl;          
            }
            break;
          } // End loop ecal associated to iTrack
        } // end case: at least one ecal element associated to iTrack
          
    } // end loop on tracks associated to hcal element iHcal

    
    // test compatibility between calo and tracker. //////////////

    double caloEnergy = -1.;
    double slopeEcal = 1.0;
    double calibEcal = std::max(0.,totalEcal);
    double calibHcal = std::max(0.,totalHcal);
    if ( newCalib_ ) {
      // caloEnergy =
        //clusterCalibration_->
        //getCalibratedEnergy(totalEcal, totalHcal, 
        //                  hclusterref->positionREP().Eta(),
        //                  hclusterref->positionREP().Phi());
      clusterCalibration_->
        getCalibratedEnergyEmbedAInHcal(calibEcal, calibHcal,
                                        hclusterref->positionREP().Eta(),
                                        hclusterref->positionREP().Phi());
      caloEnergy = calibEcal+calibHcal;
      if ( totalEcal > 0.) slopeEcal = calibEcal/totalEcal;
    } else { 
      if( totalEcal>0) { 
        caloEnergy = calibration_->energyEmHad( totalEcal, totalHcal );
        slopeEcal = calibration_->paramECALplusHCAL_slopeECAL();
        calibEcal = totalEcal * slopeEcal;
      } else { 
        caloEnergy = calibration_->energyHad( totalHcal );
        calibEcal = totalEcal;
      }
    }

    assert(caloEnergy>=0);

    //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 );
    Caloresolution *= totalChargedMomentum;
    // that of the charged particles linked to the cluster!
    double 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_;



    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
        
      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 -> loop and compute mva"<<endl;
      }
        
      // IE maxMvaIe = sortedTracks.end(); 
      reco::PFClusterRef maxMvaEcalRef;
      reco::PFClusterRef maxPSClusterRef;
      double maxMva = -PFCandidate::bigMva_;
      bool PSsignature = false;
      unsigned maxMvaiEcal= 9999;       
      unsigned maxPSiEcal= 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() ); 
          
        std::multimap< unsigned, std::pair<double, unsigned> >::iterator 
          iae = associatedEcals.find(iTrack);
          
        if( iae == associatedEcals.end() ) continue;
          
        double chi2ECAL = 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() );
          
        chi2ECAL_ = static_cast<float>( chi2ECAL );
        ptTrack_ = static_cast<float>( trackRef->pt() );
          
        float pTrack = static_cast<float>( trackRef->p() );
        float eECAL = static_cast<float>( clusterRef->energy() );
        eECALOverpTrack_ = eECAL / pTrack;
          
        double value = -PFAlgo::maxMvaCut_;
        if(mvaCut_ < PFAlgo::maxMvaCut_ ) {
          value = mergedPhotonsMVA_->EvaluateMVA( "MVA" );
          if( value>maxMva ) {
            maxMva = value;
            maxMvaEcalRef = clusterRef;
            maxMvaiEcal = iEcal;
          }
        }         
          
        if(debug_) {
          cout<<"\t\t"<<elements[iTrack]<<endl;
          cout<<"\t\t\tassociated ECAL "<<elements[iEcal]<<endl;
          cout<<"\t\t\tmva evaluation : "<<endl;
          cout<<"\t\t\t\teECALOverpTrack = "<<eECALOverpTrack_<<endl;
          cout<<"\t\t\t\tptTrack         = "<<ptTrack_<<endl;
          cout<<"\t\t\t\tchi2ECAL        = "<<chi2ECAL_<<endl;
          cout<<"\t\t\t ==>        mva = "<<value<<endl;
            
        }
        // look for PS1 elements associated to iEcal
        bool PS1found = false;
        double PS1energy = -1.;
        std::multimap<double, unsigned> sortedPS1;
        block.associatedElements( iEcal,  linkData,
                                  sortedPS1,
                                  reco::PFBlockElement::PS1 );
          
        if(debug_) cout<<"\t\t\tnumber of PS1 elements linked to this ecal: "
                       <<sortedPS1.size()<<endl;
          
        if( !sortedPS1.empty() ) 
          { // start case: at least one PS1 element associated to iEcal
            //double chi2PS1 = sortedPS1.begin()->first; //Alex//not used
            unsigned iPS1  = sortedPS1.begin()->second;
            if(debug_) cout<<"\t\t\tclosest: "
                           <<elements[iPS1]<<endl;
              
            PFBlockElement::Type type = elements[ iPS1 ].type();
            assert( type == PFBlockElement::PS1 );
            PFClusterRef PSclusterref = elements[iPS1].clusterRef();
            assert(!PSclusterref.isNull() );
              
            PS1energy = PSclusterref->energy();
            if (PS1energy > PSCut_) PS1found = true;
          } // end case: at least one PS1 element associated to iEcal
          
        // look for PS2 elements associated to iEcal
        bool PS2found = false;
        double PS2energy = -1.;
        std::multimap<double, unsigned> sortedPS2;
        block.associatedElements( iEcal,  linkData,
                                  sortedPS2,
                                  reco::PFBlockElement::PS2 );
          
        if(debug_) cout<<"\t\t\tnumber of PS2 elements linked to this ecal: "
                       <<sortedPS2.size()<<endl;
          
        if( !sortedPS2.empty() ) 
          { // start case: at least one PS2 element associated to iEcal
            //double chi2PS2 = sortedPS2.begin()->first; //Alex//noused
            unsigned iPS2  = sortedPS2.begin()->second;
            if(debug_) cout<<"\t\t\tclosest: "
                           <<elements[iPS2]<<endl;
              
            PFBlockElement::Type type = elements[ iPS2 ].type();
            assert( type == PFBlockElement::PS2 );
            PFClusterRef PSclusterref = elements[iPS2].clusterRef();
            assert(!PSclusterref.isNull() );
              
            PS2energy = PSclusterref->energy();
            if (PS2energy > PSCut_) PS2found = true;
              
          } // end case: at least one PS2 element associated to iEcal
          
          
        // evaluate PS signature
        if (PSsignature) continue;
        //PSsignature = PS1found&&PS2found;
        PSsignature = PS1found||PS2found;
        maxPSClusterRef = clusterRef;
        maxPSiEcal = iEcal;
        if (debug_) {
          cout<<" PS signature "<<PSsignature 
              <<" PS1 energy "<<PS1energy
              <<" PS2 energy "<<PS2energy<<endl;
        }
          
      }  // 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;
      vector<unsigned> elementIndices;

      // double particleEnergy = -1;
      // float  ecalEnergy = 0;
      // float  hcalEnergy = 0;
      
      if( !maxPSClusterRef.isNull() && 
          // PJ This is probably a bug too.
          // We need to save a neutral hadron if 
          // ePhoton is larger than EEcal!
          PSsignature) { //  overwrites MVA decision
        if(debug_) {
          cout<<" PS signature is true"
              <<", reconstruct a photon"<<endl;
        }
        particleEnergy.push_back(ePhoton);
        ecalEnergy.push_back(ePhoton);
        hcalEnergy.push_back(0.);
        pivotalClusterRef.push_back(maxPSClusterRef);
        iPivotal.push_back(maxPSiEcal);
      } // end case: PS signature is true
      else { // start case PS signature is false
        if( !maxMvaEcalRef.isNull() && 
            maxMva > mvaCut_) { //  start case: maxMva > mvaCut_

          if ( ePhoton < totalEcal || eNeutralHadron-calibEcal < 1E-10 ) { // PJ !!! waiting for a better MVA ;-)
            if(debug_) {
              cout<<"\t\tmaxMva = "<<maxMva<<">"<< mvaCut_
                  <<", reconstruct a photon"<<endl;
              
            }
            
            // PJ BUG !
            //if( ePhoton>totalEcal ) 
            //  ePhoton = totalEcal;
            
            pivotalClusterRef.push_back(maxMvaEcalRef);
            particleEnergy.push_back(ePhoton);
            ecalEnergy.push_back(ePhoton);
            hcalEnergy.push_back(0.);
            iPivotal.push_back(maxMvaiEcal);
            
          } else { 
            
            // Assign the ECAL energy to the photon
            pivotalClusterRef.push_back(maxMvaEcalRef);
            particleEnergy.push_back(totalEcal);
            ecalEnergy.push_back(totalEcal);
            hcalEnergy.push_back(0.);
            iPivotal.push_back(maxMvaiEcal);

            // Assign the remaining excess to a neutral hadron
            pivotalClusterRef.push_back(hclusterref);
            particleEnergy.push_back(eNeutralHadron-calibEcal);
            ecalEnergy.push_back(0.);
            hcalEnergy.push_back(eNeutralHadron-calibEcal);
            iPivotal.push_back(iHcal);

            mergedNeutralHadronEnergy = eNeutralHadron-calibEcal;

          }

        }   // end case: maxMva > mvaCut_
        else {   //  start case : maxMva < mvaCut_
          if(debug_) {
            cout<<"\t\treconstruct a hadron"<<endl;
          }
          
          // if ( totalEcal > 0. ) { 
          if ( ePhoton < totalEcal || eNeutralHadron-calibEcal <= 1E-10 ) { // PJ !!! waiting for a better MVA ;-)
            // Assign the ECAL energy to the photon
            pivotalClusterRef.push_back(maxMvaEcalRef);
            particleEnergy.push_back(ePhoton);
            ecalEnergy.push_back(ePhoton);
            hcalEnergy.push_back(0.);
            iPivotal.push_back(maxMvaiEcal);
            
          } else { 
            
            // Assign the ECAL energy to the photon
            if ( totalEcal > 0. ) { 
              pivotalClusterRef.push_back(maxMvaEcalRef);
              particleEnergy.push_back(totalEcal);
              ecalEnergy.push_back(totalEcal);
              hcalEnergy.push_back(0.);
              iPivotal.push_back(maxMvaiEcal);
            }

            // Assign the remaining excess to a neutral hadron
            pivotalClusterRef.push_back(hclusterref);
            particleEnergy.push_back(eNeutralHadron-calibEcal);
            ecalEnergy.push_back(0.);
            hcalEnergy.push_back(eNeutralHadron-calibEcal);
            iPivotal.push_back(iHcal);

            /*
              pivotalClusterRef.push_back(hclusterref);
              particleEnergy.push_back(eNeutralHadron);
              ecalEnergy.push_back(0.);
              hcalEnergy.push_back(eNeutralHadron);
              iPivotal.push_back(iHcal);
            */
            
            // keep track globally of the merged neutral hadron
            mergedNeutralHadronEnergy = eNeutralHadron-calibEcal;
          } 
        }   // end case : maxMva < mvaCut_
      }   // end case PS signature is false 


      // 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;

        unsigned tmpi = reconstructCluster( *pivotalClusterRef[iPivot], 
                                            particleEnergy[iPivot] ); 
      
        (*pfCandidates_)[tmpi].setEcalEnergy( ecalEnergy[iPivot] );
        (*pfCandidates_)[tmpi].setHcalEnergy( hcalEnergy[iPivot] );
        (*pfCandidates_)[tmpi].setPs1Energy( -1 );
        (*pfCandidates_)[tmpi].setPs2Energy( -1 );
        (*pfCandidates_)[tmpi].set_mva_nothing_gamma( maxMva );
        //       (*pfCandidates_)[tmpi].addElement(&elements[iPivotal]);
        (*pfCandidates_)[tmpi].addElementInBlock(blockref, iPivotal[iPivot]);
        // add PS elements to be done

      }

    } // excess of energy

    //COLIN will modify PFCandidate hcal (and ecal) energy
    else if(clusterRecovery_) { //  caloEnergy < totalChargedMomentum 
      // case 3: caloEnergy < totalChargedMomentum - nsigma*TotalError
      // energy is missing : look for additional hcal clusters?
      if(debug_) {
        cout<<"\t\tcase 3: CLUSTER RECOVERY "
            <<caloEnergy<<" < "
            <<totalChargedMomentum<<" - "<<nsigma
            <<" x "<<TotalError<<endl;
      }
                  
      bool investigateOtherBlocks = false;
          
      //More energy needed, looking for hcal clusters in the vicinity
      //of the tracks extrapolated to hcal in the present block
          
      if( recoveringClusters( totalChargedMomentum, 
                              totalEcal,
                              totalHcal, 
                              associatedTracks,
                              blockref,
                              active ) == KEEP_SEARCHING ){
            
        if( debug_ ){
          cout<<"\t\t\tNo satellite cluster found in current block."<<endl;
        }
        investigateOtherBlocks = true;
      }
          
      //More energy needed, looking for hcal clusters in the vicinity
      //of the tracks extrapolated to hcal in other blocks
      if( investigateOtherBlocks ){
            
        if( debug_ )
          cout<<"\t\t\tInvestigating other blocks now "<<endl;
            
        //loop on blocks
        for( IBR ih = hcalBlockRefs.begin(); 
             ih!=hcalBlockRefs.end();) {
                      
          vector<bool> active;
          active.push_back(true);
          RecoveryStatus recovery =
            recoveringClusters( totalChargedMomentum, 
                                totalEcal,
                                totalHcal, 
                                associatedTracks,
                                *ih,
                                active );
              
          if( recovery == DONE ){
            if( debug_ )
              cout<<"Cluster Recovery is done, stop here"<<endl;
            break;
          }
              
          if( !active[0] ) {
            // cluster has been used and locked in the cluster recovery
            // removing the corresponding sincle cluster block 
            // from the list of block refs to be used 
            // in cluster recovery
            // erase automatically increments the iterator
            ih = hcalBlockRefs.erase( ih );
          }
          else ++ih; // incrementing manually to continue the loop
              

        } //loop on blocks   
      } //investigating blocks
      
      // recompute calo energy after cluster recovery

      double slopeEcal = 1.0;
      double calibEcal = std::max(0.,totalEcal);
      double calibHcal = std::max(0.,totalHcal);
      if ( newCalib_ ) {
        clusterCalibration_->
          getCalibratedEnergyEmbedAInHcal(calibEcal, calibHcal,
                                          hclusterref->positionREP().Eta(),
                                          hclusterref->positionREP().Phi());
        caloEnergy = calibEcal+calibHcal;
        if ( totalEcal > 0. ) slopeEcal = calibEcal/totalEcal;
      } else { 
        if( totalEcal>0) { 
          caloEnergy = calibration_->energyEmHad( totalEcal, totalHcal );
          slopeEcal = calibration_->paramECALplusHCAL_slopeECAL();
          calibEcal = totalEcal * slopeEcal;
          calibHcal = caloEnergy - calibEcal;
        } else { 
          caloEnergy = calibration_->energyHad( totalHcal );
          calibEcal = totalEcal;
          calibHcal = caloEnergy - calibEcal;
        }
      }

    } // end if cluster recovery

    // will now share the hcal energy between the various charged hadron 
    // candidates, taking into account the potential neutral hadrons
    
    double totalHcalEnergyCalibrated = calibHcal;
    /*
    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;

    // 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;
      chargedHadron.setHcalEnergy( fraction * totalHcalEnergyCalibrated );
      // Alex: add recovered Hcal and Ecal elements?
      // create a vector <unsigned> hcalrecov
      // create a vector <unsigned> ecalrecov
      //for (unsigned i; i<hcalrecov.size(); i++){
      //unsigned iHcal = hcalrecov{i];
      //chargedHadron.addElement(&elements[iHcal]);
      //}
      //for (unsigned i; i<ecalrecov.size(); i++){
      //unsigned iEcal = ecalrecov{i];
      //chargedHadron.addElement(&elements[iEcal]);
      //}

          
    }

  } // 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;
  }
  

  for(unsigned ihcluster=0; ihcluster<hcalIs.size(); ihcluster++) {
    unsigned iHcal = hcalIs[ihcluster];
    
    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.;
    for(IE ie = ecalElems.begin(); ie != ecalElems.end(); ++ie ) {
      
      unsigned iEcal = ie->second;
      double dist = ie->first;
      double chi2 = block.chi2(iHcal,iEcal,linkData,reco::PFBlock::LINKTEST_ALL);
      PFBlockElement::Type type = elements[iEcal].type();
      assert( type == PFBlockElement::ECAL );
      
      // Check if already used
      if( !active[iEcal] ) continue;

      // Check the chi2
      if ( chi2 > 10. ) continue;

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

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

      totalEcal += eclusterRef->energy();

      active[iEcal] = false;

    }// End loop ECAL

    
    PFClusterRef hclusterRef 
      = elements[iHcal].clusterRef();
    assert( !hclusterRef.isNull() );
    
    // HCAL energy
    double totalHcal =  hclusterRef->energy();
    // Calibration
    double caloEnergy = totalHcal;
    double slopeEcal = 1.0;
    double calibEcal = totalEcal > 0. ? totalEcal : 0.;
    double calibHcal = std::max(0.,totalHcal);
    if ( newCalib_ ) {
      clusterCalibration_->
        getCalibratedEnergyEmbedAInHcal(calibEcal, calibHcal,
                                        hclusterRef->positionREP().Eta(),
                                        hclusterRef->positionREP().Phi());
      if ( calibEcal == 0. ) calibEcal = totalEcal;
      if ( calibHcal == 0. ) calibHcal = totalHcal;
      caloEnergy = calibEcal+calibHcal;
      if ( totalEcal > 0. ) slopeEcal = calibEcal/totalEcal;
    } else { 
      if( totalEcal>0) { 
        caloEnergy = calibration_->energyEmHad( totalEcal, totalHcal );
        slopeEcal = calibration_->paramECALplusHCAL_slopeECAL();
        calibEcal = totalEcal * slopeEcal;
        calibHcal = caloEnergy - calibEcal;
      } else if  (  hclusterRef->layer() != PFLayer::HCAL_HF ) { 
        caloEnergy = calibration_->energyHad( totalHcal );
        calibEcal = totalEcal;
        calibHcal = caloEnergy;
      }
    }

   // ECAL energy : calibration

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

    unsigned tmpi = reconstructCluster( *hclusterRef, 
                                        particleEnergy ); 

    
    (*pfCandidates_)[tmpi].setEcalEnergy( totalEcal );
    (*pfCandidates_)[tmpi].setHcalEnergy( calibHcal );
    (*pfCandidates_)[tmpi].setPs1Energy( -1 );
    (*pfCandidates_)[tmpi].setPs2Energy( -1 );
    (*pfCandidates_)[tmpi].addElementInBlock( blockref, iHcal );
      
  }//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;

    float ecalEnergy = calibration_->energyEm( clusterref->energy() );
    double particleEnergy = ecalEnergy;
    
    unsigned tmpi = reconstructCluster( *clusterref, 
                                        particleEnergy );
 
    (*pfCandidates_)[tmpi].setEcalEnergy( ecalEnergy );
    (*pfCandidates_)[tmpi].setHcalEnergy( 0 );
    (*pfCandidates_)[tmpi].setPs1Energy( -1 );
    (*pfCandidates_)[tmpi].setPs2Energy( -1 );
    (*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;

  double px = track.px();
  double py = track.py();
  double pz = track.pz();
  double energy = track.p();

  if( muonRef.isNonnull() ) {
    
    px = muonRef->px();
    py = muonRef->py();
    pz = muonRef->pz();
    energy = muonRef->energy();  

  }

  math::XYZTLorentzVector momentum(px,py,pz,energy);

  reco::PFCandidate::ParticleType particleType 
    = reco::PFCandidate::h;

  pfCandidates_->push_back( PFCandidate( charge, 
                                         momentum,
                                         particleType ) );
  
  pfCandidates_->back().setVertex( track.vertex() );
  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 (muonRef->isGlobalMuon() ) {
      particleType = reco::PFCandidate::mu;
      pfCandidates_->back().setParticleType( particleType );
      if (debug_) cout << "PFAlgo: particle type set to muon" << endl; 
    }
  }

  // nuclear
  if( particleType != reco::PFCandidate::mu )
    if( eltTrack->trackType(reco::PFBlockElement::T_FROM_NUCL)) {
      pfCandidates_->back().setFlag( reco::PFCandidate::T_FROM_NUCLINT, true);
      pfCandidates_->back().setNuclearRef( eltTrack->nuclearRef() );
    }
    else if( eltTrack->trackType(reco::PFBlockElement::T_TO_NUCL)) {
      pfCandidates_->back().setFlag( reco::PFCandidate::T_TO_NUCLINT, true);
      pfCandidates_->back().setNuclearRef( eltTrack->nuclearRef() );
    }

  // 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) {
  
  reco::PFCandidate::ParticleType particleType = reco::PFCandidate::X;

  // need to convert the math::XYZPoint data member of the PFCluster class=
  // to a displacement vector: 
  ROOT::Math::DisplacementVector3D< ROOT::Math::Cartesian3D<double> , ROOT::Math::DefaultCoordinateSystemTag > clusterPos( cluster.position().X(), cluster.position().Y(),cluster.position().Z() );
  

  clusterPos = clusterPos.Unit();
  clusterPos *= particleEnergy;

  // clusterPos is now a vector along the cluster direction, 
  // with a magnitude equal to the cluster energy.
  
  double mass = 0;
  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::HCAL_HF:
    particleType = PFCandidate::h_HF;
    break;
  default:
    assert(0);
  }


  ROOT::Math::LorentzVector<ROOT::Math::PxPyPzM4D<double> > 
    momentum( clusterPos.X(), clusterPos.Y(), clusterPos.Z(), mass); 

  int charge = 0;

  // mathcore is a piece of #$%
  math::XYZTLorentzVector  tmp;

  // implicit constructor not allowed
  tmp = momentum;

  pfCandidates_->push_back( PFCandidate( charge, 
                                         tmp, 
                                         particleType ) );

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

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

}




double PFAlgo::neutralHadronEnergyResolution( double clusterEnergyHCAL ) 
  const {


  return  newCalib_ ? 
    1.40/sqrt(clusterEnergyHCAL) +5.00/clusterEnergyHCAL:
    1.50/sqrt(clusterEnergyHCAL) +3.00/clusterEnergyHCAL;
  // 1.50/sqrt(clusterEnergyHCAL) +0.00/clusterEnergyHCAL;
    // 5 GeV
    // + 2.34/clusterEnergyHCAL;
    // 9 GeV
    //+ 3.15/clusterEnergyHCAL;
    // 16 GeV
    // + 4.20/clusterEnergyHCAL;
  // PJ Portnainwak!
  //  + 6.62527e-03*sqrt(clusterEnergyHCAL) -6.33966e-02; 
  // return  1.2/sqrt(clusterEnergyHCAL) + 0.05;
}




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<<"mvaCut_         "<<algo.mvaCut_<<endl;
  out<<"PSCut_          "<<algo.PSCut_<<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;
}

PFAlgo::EnergyTest 
PFAlgo::energyTest( double energyTrk, 
                    double energyHcal, 
                    double energyEcal,
                    double eta, 
                    double phi,
                    double& neutralHadron ) {
  
  //This function can be used to Test whether the calorimetric
  //energy deposition (after calibration) is compatible with 
  //the track(s) momemtum. 
  //Compatibility is defined as having the calorimetric energy
  //within +/- 3*sigma(HCAL) around the track momentum 

  double slopeEcal = 1.0;
  double calibEcal = std::max(0.,energyEcal);
  double calibHcal = std::max(0.,energyHcal);
  double eCalib = -1;
  if ( newCalib_ ) {
    clusterCalibration_->
      getCalibratedEnergyEmbedAInHcal(calibEcal, calibHcal, eta, phi);
    eCalib = calibEcal+calibHcal;
    if ( energyEcal > 0. ) slopeEcal = calibEcal/energyEcal;
  } else { 
    if( energyEcal>0) { 
      eCalib = calibration_->energyEmHad( energyEcal, energyHcal );
      slopeEcal = calibration_->paramECALplusHCAL_slopeECAL();
      calibEcal = energyEcal * slopeEcal;
      calibHcal = eCalib - calibEcal;
    } else { 
      eCalib = calibration_->energyHad( energyHcal );
      calibEcal = energyEcal;
      calibHcal = eCalib - calibEcal;
    }
  }

  /*
  double eCalib = -1;
  if( energyEcal>0)
    eCalib = calibration_->energyEmHad( energyEcal, energyHcal );
  else 
    eCalib = calibration_->energyHad( energyHcal );
  */
  
  assert(eCalib>=0);

  // PJ The resolution is on the expected charged energy
  // double reso   = neutralHadronEnergyResolution( eCalib );
  // reso *= eCalib;
  double reso   = neutralHadronEnergyResolution( energyTrk );
  reso *= energyTrk;
  
  // Alex: one should add TK errors as in main PFAlgo
  // double TotalError = sqrt(sumpError2 + Caloresolution*Caloresolution);
  //PFAlgo::recoveringClusters( double  totalChargedMomentum, double sumpError2,
  //PFAlgo::energyTest( double energyTrk, double sumpError2,  

  
  double lowerlimit = energyTrk  - nSigmaHCAL_*reso;
  double upperlimit = energyTrk  + nSigmaHCAL_*reso;

  if( debug_ ) {
    cout<<"\t\t\tEnergy test: " 
        <<"Ecal ="<<energyEcal<<" Hcal="<<energyHcal
        <<" ChargedEnergy= "<<energyTrk 
        <<" Ecalib="<<eCalib 
        <<" Resolution="<<reso 
        <<" Interval=" 
        <<lowerlimit<<" " 
        <<upperlimit<<endl;
  }

  if( (eCalib >= lowerlimit) && 
      (eCalib <= upperlimit) )
    return COMPATIBLE; //compatible
  
  if( eCalib < lowerlimit )
    return NOT_ENOUGH_ENERGY; //more energy needed

  if( eCalib > upperlimit ) {
    neutralHadron = eCalib - energyTrk - calibEcal + energyEcal;
    return TOO_MUCH_ENERGY; // too much calo energy already!
  }

  return COMPATIBLE;
}// energy Test


bool PFAlgo::isSatelliteCluster( const reco::PFRecTrack& track, 
                                 const PFCluster& cluster ){
  //This is to check whether a HCAL cluster could be 
  //a satellite cluster of a given track (charged hadron)
  //Definitions:
  //Track-Hcal: An Hcal cluster can be associated to a 
  //            track if it is found in a region around (dr<0.17 ~ 3x3 matrix)
  //            the position of the track extrapolated to the Hcal.

  bool link = false;

  if( cluster.layer() == PFLayer::HCAL_BARREL1 ||
      cluster.layer() == PFLayer::HCAL_ENDCAP ) { //Hcal case
    
    const reco::PFTrajectoryPoint& atHCAL 
      = track.extrapolatedPoint( reco::PFTrajectoryPoint::HCALEntrance );
              
    if( atHCAL.isValid() ){ //valid extrapolation?
      double tracketa = atHCAL.position().Eta();
      double trackphi = atHCAL.position().Phi();
      double hcaleta  = cluster.positionREP().Eta();
      double hcalphi  = cluster.positionREP().Phi();
                  
      //distance track-cluster
      double deta = hcaleta - tracketa;
      double dphi = acos(cos(hcalphi - trackphi));
      double dr   = sqrt(deta*deta + dphi*dphi);

      if( debug_ ){
        cout<<"\t\t\tSatellite Test " 
            <<tracketa<<" "<<trackphi<<" "
            <<hcaleta<<" "<<hcalphi<<" dr="<<dr 
            <<endl;
      }
      
      //looking if cluster is in the  
      //region around the track. 
      //Alex: Will have to adjust this cut?
      if( dr < 0.17 ) link = true;
    }//extrapolation
    
  }//HCAL
  
  return link; 
}

PFAlgo::RecoveryStatus 
PFAlgo::recoveringClusters( double  totalChargedMomentum, 
                            double  totalEcal,
                            double& totalHcal, 
                            const std::vector< reco::PFRecTrackRef >& associatedTracks,
                            const reco::PFBlockRef& blockref,
                            vector<bool>& active ) {
  

  //This fonction should be used to retrieve the satellite 
  //HCAL clusters created from interaction(s) of a charged
  //particle into the ECAL.
  
  //Note: Alex 19/11/07: Should also recover missing energy
  //in ECAL => to be implemented!

  if(debug_) 
    cout<<"\t\t\trecovering clusters"<<endl;

  //Retrieving elements of the block
  const reco::PFBlock& block
    = *blockref;
  const edm::OwnVector< reco::PFBlockElement >& 
    elements_h = block.elements();


  assert( active.size() ==  elements_h.size() );

  RecoveryStatus status = KEEP_SEARCHING;

  for(unsigned iEle=0; iEle<elements_h.size(); iEle++) {

    PFBlockElement::Type type = elements_h[iEle].type();      
    if( type != PFBlockElement::HCAL ) continue;
    
    if(debug_)
      cout<<"\t\t\thcal "<<iEle<<" ";

    //is this cluster active ?    

    if( !active[iEle] ) {
      if( debug_ ) {
        cout<<" not active. skip"<<endl;
      }
      continue;
    }    
    
    if( debug_ ) {
      cout<<" active. Is it a satellite? "<<endl;
    }

    //retrieving hcal cluster 
    const reco::PFClusterRef hclusterref 
      = elements_h[ iEle ].clusterRef();
    assert( !hclusterref.isNull() );
    
    const reco::PFCluster&  hcal = *hclusterref; 

    // if cluster active, 
    // test if it is a satellite cluster
    bool tobeadded = false;    
    for(unsigned iT=0; iT < associatedTracks.size(); ++iT) {

      assert( !associatedTracks[iT].isNull() );

      const reco::PFRecTrack& track_h 
        = *associatedTracks[iT];     

      if( isSatelliteCluster( track_h,
                              *hclusterref ) ) {
        if( debug_ ) {
          cout<<"\t\t\tsatellite of track : "<<iT<<endl;
        }
        tobeadded = true;
        break;
      }
    }//loop tracks

    //not satellite cluster?
    if( !tobeadded ) continue;
    
    totalHcal += hcal.energy();
    double neutralHadron = 0.;
    double etaH = hclusterref->positionREP().Eta();
    double phiH = hclusterref->positionREP().Phi();
    EnergyTest entest = energyTest( totalChargedMomentum,
                                    totalHcal, totalEcal,
                                    etaH, phiH,
                                    neutralHadron);
    
    if( debug_ ){
      cout<<"\t\t\tyes, adding hcal energy (" 
          <<hcal.energy()<<") to totHcal= "
          <<totalHcal 
          <<" test="<<entest<<endl;
    }

    switch( entest ) {
    case COMPATIBLE:
      //if entest ==0 then there is enough energy stop here
      if( debug_ ) cout<<"\t\t\tRecoveryStop"<<endl; 
      active[iEle]=false; 
      status = DONE;
      break;
    case NOT_ENOUGH_ENERGY:
      //if entest ==1 not enough energy continue.
      if(debug_) cout<<"\t\t\tRecoveryContinue"<<endl;
      active[iEle]=false; 
      status = KEEP_SEARCHING;   
      continue;
    case TOO_MUCH_ENERGY:
      //if entest ==2 too much energy has been added, 
      //create a neutral hadron
      unsigned tmpi = reconstructCluster( *hclusterref, 
                                          neutralHadron ); 

      // std::cout << "A cluster recovered here " << neutralHadron << std::endl;
      
      (*pfCandidates_)[tmpi].setEcalEnergy( 0. );
      (*pfCandidates_)[tmpi].setHcalEnergy( neutralHadron );
      (*pfCandidates_)[tmpi].setPs1Energy( -1 );
      (*pfCandidates_)[tmpi].setPs2Energy( -1 );
      (*pfCandidates_)[tmpi].addElementInBlock( blockref, iEle );
       
      active[iEle]=false; 
      status = DONE;
      break;
    default:
      assert(0);
      //      continue;
    }

  }//loop elements

  return status;
} 

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

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

---