#include <PFAlgo.h>

Inheritance diagram for PFAlgo:

Public Member Functions
void	checkCleaning (const reco::PFRecHitCollection &cleanedHF)
	Check HF Cleaning. More...

	PFAlgo ()
	constructor More...

const std::auto_ptr < reco::PFCandidateCollection > &	pfCandidates () const

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

void	reconstructParticles (const reco::PFBlockHandle &blockHandle)

virtual void	reconstructParticles (const reco::PFBlockCollection &blocks)
	reconstruct particles More...

void	setAlgo (int algo)

void	setCandConnectorParameters (const edm::ParameterSet &iCfgCandConnector)

void	setCandConnectorParameters (bool bCorrect, bool bCalibPrimary, double dptRel_PrimaryTrack, double dptRel_MergedTrack, double ptErrorSecondary, std::vector< double > nuclCalibFactors)

void	setDebug (bool debug)

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

void	setEGElectronCollection (const reco::GsfElectronCollection &egelectrons)

void	setElectronExtraRef (const edm::OrphanHandle< reco::PFCandidateElectronExtraCollection > &extrah)

void	setHOTag (bool ho)

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

void	setPFEleParameters (double mvaEleCut, std::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=false, bool usePFSCEleCalib=true, bool useEGElectrons=false, bool useEGammaSupercluster=true)

void	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)

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

void	setPFPhotonRegWeights (const GBRForest LCorrForestEB, const GBRForest LCorrForestEE, const GBRForest GCorrForestBarrel, const GBRForest GCorrForestEndcapHr9, const GBRForest GCorrForestEndcapLr9, const GBRForest PFEcalResolution)

void	setPFVertexParameters (bool useVertex, const reco::VertexCollection &primaryVertices)

void	setPhotonExtraRef (const edm::OrphanHandle< reco::PFCandidatePhotonExtraCollection > &pf_extrah)

void	setPostHFCleaningParameters (bool postHFCleaning, double minHFCleaningPt, double minSignificance, double maxSignificance, double minSignificanceReduction, double maxDeltaPhiPt, double minDeltaMet)

boost::shared_ptr < PFEnergyCalibration >	thePFEnergyCalibration ()
	return the pointer to the calibration function More...

std::auto_ptr < reco::PFCandidateCollection >	transferAddedMuonCandidates ()

std::auto_ptr < reco::PFCandidateCollection >	transferCandidates ()

std::auto_ptr < reco::PFCandidateCollection > &	transferCleanedCandidates ()

std::auto_ptr < reco::PFCandidateCollection >	transferCleanedTrackerAndGlobalMuonCandidates ()

std::auto_ptr < reco::PFCandidateCollection >	transferCosmicsMuonCleanedCandidates ()

std::auto_ptr < reco::PFCandidateCollection >	transferElectronCandidates ()

std::auto_ptr < reco::PFCandidateElectronExtraCollection >	transferElectronExtra ()

std::auto_ptr < reco::PFCandidateCollection >	transferFakeMuonCleanedCandidates ()

std::auto_ptr < reco::PFCandidatePhotonExtraCollection >	transferPhotonExtra ()

std::auto_ptr < reco::PFCandidateCollection >	transferPunchThroughHadronCleanedCandidates ()

std::auto_ptr < reco::PFCandidateCollection >	transferPunchThroughMuonCleanedCandidates ()

virtual	~PFAlgo ()
	destructor More...

Protected Member Functions
void	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)
	Associate PS clusters to a given ECAL cluster, and return their energy. More...

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

double	neutralHadronEnergyResolution (double clusterEnergy, double clusterEta) const
	todo: use PFClusterTools for this More...

double	nSigmaHCAL (double clusterEnergy, double clusterEta) const

void	postCleaning ()

virtual void	processBlock (const reco::PFBlockRef &blockref, std::list< reco::PFBlockRef > &hcalBlockRefs, std::list< reco::PFBlockRef > &ecalBlockRefs)

unsigned	reconstructCluster (const reco::PFCluster &cluster, double particleEnergy, bool useDirection=false, double particleX=0., double particleY=0., double particleZ=0.)

unsigned	reconstructTrack (const reco::PFBlockElement &elt)

Protected Attributes
std::auto_ptr < reco::PFCandidateCollection >	pfAddedMuonCandidates_
	the collection of added muon candidates More...

std::auto_ptr < reco::PFCandidateCollection >	pfCandidates_

std::auto_ptr < reco::PFCandidateCollection >	pfCleanedCandidates_

std::auto_ptr < reco::PFCandidateCollection >	pfCleanedTrackerAndGlobalMuonCandidates_
	the collection of tracker/global cleaned muon candidates More...

std::auto_ptr < reco::PFCandidateCollection >	pfCosmicsMuonCleanedCandidates_
	the collection of cosmics cleaned muon candidates More...

std::auto_ptr < reco::PFCandidateCollection >	pfElectronCandidates_
	the unfiltered electron collection More...

reco::PFCandidateElectronExtraCollection	pfElectronExtra_
	the unfiltered electron collection More...

std::auto_ptr < reco::PFCandidateCollection >	pfFakeMuonCleanedCandidates_
	the collection of fake cleaned muon candidates More...

std::auto_ptr < reco::PFCandidateCollection >	pfPhotonCandidates_
	the unfiltered photon collection More...

reco::PFCandidatePhotonExtraCollection	pfPhotonExtra_
	the extra photon collection More...

std::auto_ptr < reco::PFCandidateCollection >	pfPunchThroughHadronCleanedCandidates_
	the collection of punch-through cleaned neutral hadron candidates More...

std::auto_ptr < reco::PFCandidateCollection >	pfPunchThroughMuonCleanedCandidates_
	the collection of punch-through cleaned muon candidates More...

Private Member Functions
reco::PFBlockRef	createBlockRef (const reco::PFBlockCollection &blocks, unsigned bi)

Private Attributes
int	algo_

bool	applyCrackCorrectionsElectrons_

reco::PFBlockHandle	blockHandle_
	input block handle (full framework case) More...

boost::shared_ptr < PFEnergyCalibration >	calibration_

double	coneEcalIsoForEgammaSC_

double	coneTrackIsoForEgammaSC_

PFCandConnector	connector_

bool	debug_

double	dptRel_DispVtx_

std::vector< double >	factors45_

double	maxDeltaPhiPt_

double	maxSignificance_

double	minDeltaMet_

double	minHFCleaningPt_

double	minSignificance_

double	minSignificanceReduction_

std::vector< double >	muonECAL_

std::vector< double >	muonHCAL_
	Variables for muons and fakes. More...

std::vector< double >	muonHO_

double	mvaEleCut_

std::string	mvaWeightFileEleID_
	Variables for PFElectrons. More...

double	nSigmaECAL_
	number of sigma to judge energy excess in ECAL More...

double	nSigmaHCAL_
	number of sigma to judge energy excess in HCAL More...

double	nSigmaTRACK_

unsigned int	nTrackIsoForEgammaSC_

PFElectronAlgo *	pfele_

PFPhotonAlgo *	pfpho_

bool	postHFCleaning_

bool	postMuonCleaning_

reco::Vertex	primaryVertex_

double	ptError_

bool	rejectTracks_Bad_

bool	rejectTracks_Step45_

std::vector< double >	setchi2Values_

double	sumEtEcalIsoForEgammaSC_barrel_

double	sumEtEcalIsoForEgammaSC_endcap_

double	sumPtTrackIsoForEgammaSC_barrel_

double	sumPtTrackIsoForEgammaSC_endcap_

boost::shared_ptr < PFEnergyCalibrationHF >	thepfEnergyCalibrationHF_

boost::shared_ptr < PFSCEnergyCalibration >	thePFSCEnergyCalibration_

double	useBestMuonTrack_

bool	useEGammaSupercluster_

bool	useEGElectrons_

bool	useHO_

bool	usePFConversions_

bool	usePFDecays_

bool	usePFElectrons_

bool	usePFMuonMomAssign_

bool	usePFNuclearInteractions_

bool	usePFPhotons_

bool	usePFSCEleCalib_

bool	useVertices_

Friends
std::ostream &	operator<< (std::ostream &out, const PFAlgo &algo)

Detailed Description

Definition at line 48 of file PFAlgo.h.

Constructor & Destructor Documentation

PFAlgo::PFAlgo ( )

constructor

Definition at line 57 of file PFAlgo.cc.

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

PFAlgo::~PFAlgo ( )

virtual

destructor

Definition at line 68 of file PFAlgo.cc.

References pfele_, pfpho_, usePFElectrons_, and usePFPhotons_.

                 {
   if (usePFElectrons_) delete pfele_;
   if (usePFPhotons_)     delete pfpho_;
 }

Member Function Documentation

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
	)

protected

Associate PS clusters to a given ECAL cluster, and return their energy.

Definition at line 3306 of file PFAlgo.cc.

References reco::PFBlock::associatedElements(), reco::PFBlockElement::ECAL, edm::Ref< C, T, F >::isNull(), and reco::PFBlock::LINKTEST_ALL.

Referenced by processBlock().

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

void PFAlgo::checkCleaning ( const reco::PFRecHitCollection & cleanedHF )

Check HF Cleaning.

Definition at line 3504 of file PFAlgo.cc.

References gather_cfg::cout, debug_, reco::PFRecHit::energy(), i, j, reco::PFRecHit::layer(), minDeltaMet_, GetRecoTauVFromDQM_MC_cff::next, pfCandidates_, reco::PFRecHit::position(), reco::LeafCandidate::pt(), reco::LeafCandidate::px(), reco::LeafCandidate::py(), reconstructCluster(), createPayload::skip, mathSSE::sqrt(), and reco::PFRecHit::time().

Referenced by PFRootEventManager::particleFlow().

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

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

private

create a reference to a block, transient or persistent depending on the needs

Definition at line 3294 of file PFAlgo.cc.

References blockHandle_, and edm::HandleBase::isValid().

Referenced by reconstructParticles().

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

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

protected

Definition at line 3361 of file PFAlgo.cc.

References reco::PFBlockElement::T_FROM_DISP, reco::PFBlockElement::T_FROM_V0, reco::PFBlockElement::T_TO_DISP, reco::PFBlockElement::trackType(), usePFDecays_, and usePFNuclearInteractions_.

Referenced by processBlock(), and reconstructTrack().

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

double PFAlgo::neutralHadronEnergyResolution	(	double	clusterEnergy,
		double	clusterEta
	)		const

protected

todo: use PFClusterTools for this

Returns: calibrated energy of a photon; calibrated energy of a neutral hadron, which can leave some energy in the ECAL ( energyECAL>0 )

Definition at line 3241 of file PFAlgo.cc.

References mathSSE::sqrt().

Referenced by processBlock().

                                                                                 {
 
   // 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	clusterEnergy,
		double	clusterEta
	)		const

protected

Definition at line 3256 of file PFAlgo.cc.

References create_public_lumi_plots::exp, and nSigmaHCAL_.

Referenced by processBlock(), and setParameters().

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

const std::auto_ptr< reco::PFCandidateCollection >& PFAlgo::pfCandidates ( ) const

inline

Returns: collection of candidates

Definition at line 176 of file PFAlgo.h.

References pfCandidates_.

Referenced by operator<<().

                                                                        {
     return pfCandidates_;
   }

void PFAlgo::postCleaning ( )

protected

Definition at line 3397 of file PFAlgo.cc.

References gather_cfg::cout, SiPixelRawToDigiRegional_cfi::deltaPhi, reco::PFCandidate::egamma_HF, reco::PFCandidate::h_HF, i, j, maxDeltaPhiPt_, maxSignificance_, minDeltaMet_, minHFCleaningPt_, minSignificance_, minSignificanceReduction_, GetRecoTauVFromDQM_MC_cff::next, reco::PFCandidate::particleId(), pfCandidates_, pfCleanedCandidates_, postHFCleaning_, reco::LeafCandidate::pt(), reco::LeafCandidate::px(), reco::LeafCandidate::py(), createPayload::skip, and mathSSE::sqrt().

Referenced by reconstructParticles().

                      { 
   
   // 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::postMuonCleaning	(	const edm::Handle< reco::MuonCollection > &	muonh,
		const reco::VertexCollection &	primaryVertices
	)

Definition at line 3609 of file PFAlgo.cc.

References DeDxDiscriminatorTools::charge(), gather_cfg::cout, dPhi(), PFRecoTauDiscriminationAgainstElectronDeadECAL_cfi::dR, reco::PFCandidate::elementsInBlocks(), relval_parameters_module::energy, reco::LeafCandidate::eta(), reco::PFCandidate::h, reco::PFCandidate::h0, i, PFMuonAlgo::isIsolatedMuon(), edm::Ref< C, T, F >::isNonnull(), edm::Ref< C, T, F >::key(), max(), reco::PFCandidate::mu, reco::PFCandidate::muonRef(), reco::LeafCandidate::p(), p4, reco::PFCandidate::particleId(), pfAddedMuonCandidates_, pfCandidates_, pfCleanedTrackerAndGlobalMuonCandidates_, pfCosmicsMuonCleanedCandidates_, pfFakeMuonCleanedCandidates_, pfMuons_cff::pfMuons, pfPunchThroughHadronCleanedCandidates_, pfPunchThroughMuonCleanedCandidates_, primaryVertex_, reco::LeafCandidate::pt(), reco::LeafCandidate::px(), reco::LeafCandidate::py(), CosmicsPD_Skims::radius, reco::PFCandidate::rawEcalEnergy(), reco::PFCandidate::rawHcalEnergy(), findQualityFiles::size, mathSSE::sqrt(), reco::PFCandidate::trackRef(), useBestMuonTrack_, reco::PFCandidate::vx(), reco::PFCandidate::vy(), reco::Vertex::x(), and reco::Vertex::y().

Referenced by PFRootEventManager::particleFlow().

                                                            {
   
 
   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::processBlock	(	const reco::PFBlockRef &	blockref,
		std::list< reco::PFBlockRef > &	hcalBlockRefs,
		std::list< reco::PFBlockRef > &	ecalBlockRefs
	)

protectedvirtual

process one block. can be reimplemented in more sophisticated algorithms

Reimplemented in PFAlgoTestBenchConversions, and PFAlgoTestBenchElectrons.

Definition at line 451 of file PFAlgo.cc.

Referenced by reconstructParticles().

                                                                     { 
   
   // 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::reconstructCluster	(	const reco::PFCluster &	cluster,
		double	particleEnergy,
		bool	useDirection = `false`,
		double	particleX = `0.`,
		double	particleY = `0.`,
		double	particleZ = `0.`
	)

protected

Reconstruct a neutral particle from a cluster. If chargedEnergy is specified, the neutral particle is created only if the cluster energy is significantly larger than the chargedEnergy. In this case, the energy of the neutral particle is cluster energy - chargedEnergy

Definition at line 3129 of file PFAlgo.cc.

References DeDxDiscriminatorTools::charge(), gather_cfg::cout, debug_, PFLayer::ECAL_BARREL, PFLayer::ECAL_ENDCAP, PFLayer::HCAL_BARREL1, PFLayer::HCAL_ENDCAP, PFLayer::HF_EM, PFLayer::HF_HAD, reco::PFCluster::layer(), scaleCards::mass, pfCandidates_, reco::CaloCluster::position(), primaryVertex_, mathSSE::sqrt(), tmp, useVertices_, reco::PFCandidate::X, reco::Vertex::x(), reco::Vertex::y(), and reco::Vertex::z().

Referenced by checkCleaning(), and processBlock().

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

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

reconstruct particles (full framework case) will keep track of the block handle to build persistent references, and call reconstructParticles( const reco::PFBlockCollection& blocks )

Definition at line 305 of file PFAlgo.cc.

References blockHandle_.

Referenced by PFRootEventManager::particleFlow().

                                                                         {
 
   blockHandle_ = blockHandle;
   reconstructParticles( *blockHandle_ ); 
 }

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

virtual

reconstruct particles

Definition at line 313 of file PFAlgo.cc.

References Association::block, gather_cfg::cout, createBlockRef(), debug_, reco::PFBlockElement::ECAL, reco::PFBlock::elements(), asciidump::elements, relativeConstraints::empty, reco::PFBlockElement::HCAL, reco::PFBlockElement::HO, i, pfAddedMuonCandidates_, pfCandidates_, pfCleanedCandidates_, pfCleanedTrackerAndGlobalMuonCandidates_, pfCosmicsMuonCleanedCandidates_, pfElectronCandidates_, pfElectronExtra_, pfFakeMuonCleanedCandidates_, pfPhotonCandidates_, pfPhotonExtra_, pfPunchThroughHadronCleanedCandidates_, pfPunchThroughMuonCleanedCandidates_, postCleaning(), processBlock(), and edm::OwnVector< T, P >::size().

                                                                        {
 
   // 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();
 }

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

protected

Reconstruct a charged hadron from a track Returns the index of the newly created candidate in pfCandidates_

Definition at line 2934 of file PFAlgo.cc.

References DeDxDiscriminatorTools::charge(), reco::TrackBase::charge(), gather_cfg::cout, debug_, reco::PFBlockElementTrack::displacedVertexRef(), dptRel_DispVtx_, relval_parameters_module::energy, reco::PFCandidate::h, reco::TrackBase::hitPattern(), isFromSecInt(), PFMuonAlgo::isGlobalTightMuon(), PFMuonAlgo::isIsolatedMuon(), PFMuonAlgo::isMuon(), edm::Ref< C, T, F >::isNonnull(), PFMuonAlgo::isTrackerTightMuon(), reco::PFCandidate::mu, reco::PFBlockElementTrack::muonRef(), reco::HitPattern::numberOfValidPixelHits(), reco::HitPattern::numberOfValidTrackerHits(), reco::TrackBase::p(), pfCandidates_, reco::PFBlockElementTrack::positionAtECALEntrance(), PFMuonAlgo::printMuonProperties(), reco::TrackBase::ptError(), reco::TrackBase::px(), reco::TrackBase::py(), reco::TrackBase::pz(), mathSSE::sqrt(), reco::PFBlockElement::T_FROM_DISP, reco::PFCandidate::T_FROM_DISP, reco::PFBlockElement::T_TO_DISP, reco::PFCandidate::T_TO_DISP, reco::PFBlockElementTrack::trackRef(), useBestMuonTrack_, and usePFMuonMomAssign_.

Referenced by processBlock().

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

void PFAlgo::setAlgo ( int algo )

inline

Definition at line 59 of file PFAlgo.h.

References algo, and algo_.

Referenced by PFRootEventManager::readOptions().

59 {algo_ = algo;}

PFAlgo::algo_

int algo_

Definition: PFAlgo.h:356

algo

LimitAlgo * algo

Definition: Combine.cc:60

void PFAlgo::setCandConnectorParameters ( const edm::ParameterSet & iCfgCandConnector )

inline

Definition at line 68 of file PFAlgo.h.

References connector_, and PFCandConnector::setParameters().

Referenced by PFRootEventManager::readOptions().

                                                                              {
     connector_.setParameters(iCfgCandConnector);
   }

void PFAlgo::setCandConnectorParameters	(	bool	bCorrect,
		bool	bCalibPrimary,
		double	dptRel_PrimaryTrack,
		double	dptRel_MergedTrack,
		double	ptErrorSecondary,
		std::vector< double >	nuclCalibFactors
	)

inline

Definition at line 72 of file PFAlgo.h.

References connector_, and PFCandConnector::setParameters().

                                                      {
     connector_.setParameters(bCorrect, bCalibPrimary, dptRel_PrimaryTrack, dptRel_MergedTrack, ptErrorSecondary, nuclCalibFactors);
   }

void PFAlgo::setDebug ( bool debug )

inline

Definition at line 61 of file PFAlgo.h.

References connector_, debug, debug_, and PFCandConnector::setDebug().

Referenced by PFRootEventManager::readOptions().

61 {debug_ = debug; connector_.setDebug(debug_);}

PFAlgo::connector_

PFCandConnector connector_

Definition: PFAlgo.h:398

PFAlgo::debug_

bool debug_

Definition: PFAlgo.h:357

debug

#define debug

Definition: MEtoEDMFormat.h:34

PFCandConnector::setDebug

void setDebug(bool debug)

Definition: PFCandConnector.h:64

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

Definition at line 264 of file PFAlgo.cc.

References dptRel_DispVtx_, rejectTracks_Bad_, rejectTracks_Step45_, usePFConversions_, usePFDecays_, and usePFNuclearInteractions_.

Referenced by PFRootEventManager::readOptions().

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

void PFAlgo::setEGElectronCollection ( const reco::GsfElectronCollection & egelectrons )

Definition at line 3392 of file PFAlgo.cc.

References pfele_, PFElectronAlgo::setEGElectronCollection(), and useEGElectrons_.

Referenced by PFRootEventManager::particleFlow().

                                                                              {
   if(useEGElectrons_ && pfele_) pfele_->setEGElectronCollection(egelectrons);
 }

void PFAlgo::setElectronExtraRef ( const edm::OrphanHandle< reco::PFCandidateElectronExtraCollection > & extrah )

Definition at line 4137 of file PFAlgo.cc.

References alignCSCRings::e, pfCandidates_, pfElectronCandidates_, pfElectronExtra_, findQualityFiles::size, and usePFElectrons_.

Referenced by PFRootEventManager::particleFlow().

                                                                                                        {
   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::setHOTag ( bool ho )

inline

Definition at line 58 of file PFAlgo.h.

References useHO_.

Referenced by PFRootEventManager::readOptions().

58 { useHO_ = ho;}

PFAlgo::useHO_

bool useHO_

Definition: PFAlgo.h:355

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

Definition at line 75 of file PFAlgo.cc.

References calibration_, nSigmaECAL_, nSigmaHCAL(), nSigmaHCAL_, and thepfEnergyCalibrationHF_.

Referenced by PFRootEventManager::readOptions().

                                                                                        {
 
   nSigmaECAL_ = nSigmaECAL;
   nSigmaHCAL_ = nSigmaHCAL;
 
   calibration_ = calibration;
   thepfEnergyCalibrationHF_ = thepfEnergyCalibrationHF;
 
 }

void PFAlgo::setPFEleParameters	(	double	mvaEleCut,
		std::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 = `false`,
		bool	usePFSCEleCalib = `true`,
		bool	useEGElectrons = `false`,
		bool	useEGammaSupercluster = `true`
	)

Definition at line 90 of file PFAlgo.cc.

References applyCrackCorrectionsElectrons_, coneEcalIsoForEgammaSC_, coneTrackIsoForEgammaSC_, mvaEleCut_, mvaWeightFileEleID_, nTrackIsoForEgammaSC_, pfele_, sumEtEcalIsoForEgammaSC_barrel_, sumEtEcalIsoForEgammaSC_endcap_, sumPtTrackIsoForEgammaSC_barrel_, sumPtTrackIsoForEgammaSC_endcap_, thePFSCEnergyCalibration_, useEGammaSupercluster_, useEGElectrons_, usePFElectrons_, and usePFSCEleCalib_.

Referenced by PFRootEventManager::readOptions().

                                            {
   
   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::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
	)

Definition at line 227 of file PFAlgo.cc.

References factors45_, muonECAL_, muonHCAL_, muonHO_, nSigmaTRACK_, ptError_, useBestMuonTrack_, and usePFMuonMomAssign_.

Referenced by PFRootEventManager::readOptions().

 {
   muonHCAL_ = muonHCAL;
   muonECAL_ = muonECAL;
   muonHO_ = muonHO;
   nSigmaTRACK_ = nSigmaTRACK;
   ptError_ = ptError;
   factors45_ = factors45;
   usePFMuonMomAssign_ = usePFMuonMomAssign;
   useBestMuonTrack_ = useBestMuonTrack;
 }

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

Definition at line 152 of file PFAlgo.cc.

References alignCSCRings::e, AlCaHLTBitMon_ParallelJobs::p, pfpho_, primaryVertex_, usePFPhotons_, and useVertices_.

Referenced by PFRootEventManager::readOptions().

  {
 
   usePFPhotons_ = usePFPhotons;
 
   //for MVA pass PV if there is one in the collection otherwise pass a dummy    
   reco::Vertex dummy;  
   const reco::Vertex* pv=&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);  
   }  
    
   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 );  
   }  
    
   
   pfpho_ = new PFPhotonAlgo(mvaWeightFileConvID, 
                 mvaConvCut, 
                 useReg,
                 X0_Map,  
                 *pv,
                 thePFEnergyCalibration,
                             sumPtTrackIsoForPhoton,
                             sumPtTrackIsoSlopeForPhoton
                 );
   return;
 }

void PFAlgo::setPFPhotonRegWeights	(	const GBRForest *	LCorrForestEB,
		const GBRForest *	LCorrForestEE,
		const GBRForest *	GCorrForestBarrel,
		const GBRForest *	GCorrForestEndcapHr9,
		const GBRForest *	GCorrForestEndcapLr9,
		const GBRForest *	PFEcalResolution
	)

Definition at line 214 of file PFAlgo.cc.

References pfpho_, and PFPhotonAlgo::setGBRForest().

Referenced by PFRootEventManager::readOptions().

                     {
   
   pfpho_->setGBRForest(LCorrForestEB,LCorrForestEE,
                GCorrForestBarrel, GCorrForestEndcapHr9, 
                GCorrForestEndcapLr9, PFEcalResolution);
 }

void PFAlgo::setPFVertexParameters	(	bool	useVertex,
		const reco::VertexCollection &	primaryVertices
	)

Definition at line 282 of file PFAlgo.cc.

References i, pfpho_, primaryVertex_, PFPhotonAlgo::setnPU(), usePFPhotons_, and useVertices_.

Referenced by PFRootEventManager::particleFlow().

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

void PFAlgo::setPhotonExtraRef ( const edm::OrphanHandle< reco::PFCandidatePhotonExtraCollection > & pf_extrah )

Definition at line 4189 of file PFAlgo.cc.

References newFWLiteAna::found, pfCandidates_, pfPhotonExtra_, findQualityFiles::size, and usePFPhotons_.

                                                                                                       {
   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
       }
     }
   }      
 }

void PFAlgo::setPostHFCleaningParameters	(	bool	postHFCleaning,
		double	minHFCleaningPt,
		double	minSignificance,
		double	maxSignificance,
		double	minSignificanceReduction,
		double	maxDeltaPhiPt,
		double	minDeltaMet
	)

Definition at line 247 of file PFAlgo.cc.

References maxDeltaPhiPt_, maxSignificance_, minDeltaMet_, minHFCleaningPt_, minSignificance_, minSignificanceReduction_, and postHFCleaning_.

Referenced by PFRootEventManager::readOptions().

                                         { 
   postHFCleaning_ = postHFCleaning;
   minHFCleaningPt_ = minHFCleaningPt;
   minSignificance_ = minSignificance;
   maxSignificance_ = maxSignificance;
   minSignificanceReduction_= minSignificanceReduction;
   maxDeltaPhiPt_ = maxDeltaPhiPt;
   minDeltaMet_ = minDeltaMet;
 }

boost::shared_ptr<PFEnergyCalibration> PFAlgo::thePFEnergyCalibration ( )

inline

return the pointer to the calibration function

Definition at line 244 of file PFAlgo.h.

References calibration_.

                                                                 { 
     return calibration_;
   }

std::auto_ptr< reco::PFCandidateCollection > PFAlgo::transferAddedMuonCandidates ( )

inline

Returns: collection of added muon candidates

Definition at line 234 of file PFAlgo.h.

References pfAddedMuonCandidates_.

                                                                          {
     return pfAddedMuonCandidates_;
   }

std::auto_ptr< reco::PFCandidateCollection > PFAlgo::transferCandidates ( )

inline

Returns: auto_ptr to the collection of candidates (transfers ownership)

Definition at line 239 of file PFAlgo.h.

References PFCandConnector::connect(), connector_, and pfCandidates_.

Referenced by PFRootEventManager::particleFlow().

                                                                  {
     return connector_.connect(pfCandidates_);
   }

std::auto_ptr< reco::PFCandidateCollection >& PFAlgo::transferCleanedCandidates ( )

inline

Returns: collection of cleaned HF candidates

Definition at line 204 of file PFAlgo.h.

References pfCleanedCandidates_.

                                                                         {
     return pfCleanedCandidates_;
   }

std::auto_ptr< reco::PFCandidateCollection > PFAlgo::transferCleanedTrackerAndGlobalMuonCandidates ( )

inline

Returns: collection of tracker/global cleaned muon candidates

Definition at line 214 of file PFAlgo.h.

References pfCleanedTrackerAndGlobalMuonCandidates_.

                                                                                            {
     return pfCleanedTrackerAndGlobalMuonCandidates_;
   }

std::auto_ptr< reco::PFCandidateCollection > PFAlgo::transferCosmicsMuonCleanedCandidates ( )

inline

Returns: collection of cosmics cleaned muon candidates

Definition at line 209 of file PFAlgo.h.

References pfCosmicsMuonCleanedCandidates_.

                                                                                   {
     return pfCosmicsMuonCleanedCandidates_;
   }

std::auto_ptr< reco::PFCandidateCollection> PFAlgo::transferElectronCandidates ( )

inline

Returns: the unfiltered electron collection

Definition at line 181 of file PFAlgo.h.

References pfElectronCandidates_.

                                                                         {
     return pfElectronCandidates_;
   }

std::auto_ptr< reco::PFCandidateElectronExtraCollection> PFAlgo::transferElectronExtra ( )

inline

Returns: the unfiltered electron extra collection

Definition at line 187 of file PFAlgo.h.

References pfElectronExtra_, and query::result.

Referenced by PFRootEventManager::particleFlow().

                                                                                 {
     std::auto_ptr< reco::PFCandidateElectronExtraCollection> result(new reco::PFCandidateElectronExtraCollection);
     result->insert(result->end(),pfElectronExtra_.begin(),pfElectronExtra_.end());
     return result;
   }

std::auto_ptr< reco::PFCandidateCollection > PFAlgo::transferFakeMuonCleanedCandidates ( )

inline

Returns: collection of fake cleaned muon candidates

Definition at line 219 of file PFAlgo.h.

References pfFakeMuonCleanedCandidates_.

                                                                                {
     return pfFakeMuonCleanedCandidates_;
   }

std::auto_ptr< reco::PFCandidatePhotonExtraCollection> PFAlgo::transferPhotonExtra ( )

inline

Returns: the unfiltered photon extra collection

Definition at line 196 of file PFAlgo.h.

References pfPhotonExtra_, and query::result.

                                                                             {
     std::auto_ptr< reco::PFCandidatePhotonExtraCollection> result(new reco::PFCandidatePhotonExtraCollection);
     result->insert(result->end(),pfPhotonExtra_.begin(),pfPhotonExtra_.end());
     return result;
   }

std::auto_ptr< reco::PFCandidateCollection > PFAlgo::transferPunchThroughHadronCleanedCandidates ( )

inline

Returns: collection of punch-through cleaned neutral hadron candidates

Definition at line 229 of file PFAlgo.h.

References pfPunchThroughHadronCleanedCandidates_.

                                                                                          {
     return pfPunchThroughHadronCleanedCandidates_;
   }

std::auto_ptr< reco::PFCandidateCollection > PFAlgo::transferPunchThroughMuonCleanedCandidates ( )

inline

Returns: collection of punch-through cleaned muon candidates

Definition at line 224 of file PFAlgo.h.

References pfPunchThroughMuonCleanedCandidates_.

                                                                                        {
     return pfPunchThroughMuonCleanedCandidates_;
   }

Friends And Related Function Documentation

std::ostream& operator<<	(	std::ostream &	out,
		const PFAlgo &	algo
	)

friend

Member Data Documentation

int PFAlgo::algo_

private

Definition at line 356 of file PFAlgo.h.

Referenced by setAlgo().

bool PFAlgo::applyCrackCorrectionsElectrons_

private

Definition at line 365 of file PFAlgo.h.

Referenced by setPFEleParameters().

reco::PFBlockHandle PFAlgo::blockHandle_

private

input block handle (full framework case)

Definition at line 343 of file PFAlgo.h.

Referenced by createBlockRef(), and reconstructParticles().

boost::shared_ptr<PFEnergyCalibration> PFAlgo::calibration_

private

Definition at line 351 of file PFAlgo.h.

Referenced by operator<<(), processBlock(), setParameters(), and thePFEnergyCalibration().

double PFAlgo::coneEcalIsoForEgammaSC_

private

Definition at line 371 of file PFAlgo.h.

Referenced by setPFEleParameters().

double PFAlgo::coneTrackIsoForEgammaSC_

private

Definition at line 374 of file PFAlgo.h.

Referenced by setPFEleParameters().

PFCandConnector PFAlgo::connector_

private

A tool used for a postprocessing of displaced vertices based on reconstructed PFCandidates

Definition at line 398 of file PFAlgo.h.

Referenced by setCandConnectorParameters(), setDebug(), and transferCandidates().

bool PFAlgo::debug_

private

Definition at line 357 of file PFAlgo.h.

Referenced by checkCleaning(), processBlock(), reconstructCluster(), reconstructParticles(), reconstructTrack(), and setDebug().

double PFAlgo::dptRel_DispVtx_

private

Maximal relative uncertainty on the tracks going to or incoming from the displcaed vertex to be used in the PFAlgo

Definition at line 393 of file PFAlgo.h.

Referenced by processBlock(), reconstructTrack(), and setDisplacedVerticesParameters().

std::vector<double> PFAlgo::factors45_

private

Definition at line 406 of file PFAlgo.h.

Referenced by processBlock(), and setPFMuonAndFakeParameters().

double PFAlgo::maxDeltaPhiPt_

private

Definition at line 415 of file PFAlgo.h.

Referenced by postCleaning(), and setPostHFCleaningParameters().

double PFAlgo::maxSignificance_

private

Definition at line 413 of file PFAlgo.h.

Referenced by postCleaning(), and setPostHFCleaningParameters().

double PFAlgo::minDeltaMet_

private

Definition at line 416 of file PFAlgo.h.

Referenced by checkCleaning(), postCleaning(), and setPostHFCleaningParameters().

double PFAlgo::minHFCleaningPt_

private

Definition at line 411 of file PFAlgo.h.

Referenced by postCleaning(), and setPostHFCleaningParameters().

double PFAlgo::minSignificance_

private

Definition at line 412 of file PFAlgo.h.

Referenced by postCleaning(), and setPostHFCleaningParameters().

double PFAlgo::minSignificanceReduction_

private

Definition at line 414 of file PFAlgo.h.

Referenced by postCleaning(), and setPostHFCleaningParameters().

std::vector<double> PFAlgo::muonECAL_

private

Definition at line 402 of file PFAlgo.h.

Referenced by processBlock(), and setPFMuonAndFakeParameters().

std::vector<double> PFAlgo::muonHCAL_

private

Variables for muons and fakes.

Definition at line 401 of file PFAlgo.h.

Referenced by processBlock(), and setPFMuonAndFakeParameters().

std::vector<double> PFAlgo::muonHO_

private

Definition at line 403 of file PFAlgo.h.

Referenced by processBlock(), and setPFMuonAndFakeParameters().

double PFAlgo::mvaEleCut_

private

Definition at line 362 of file PFAlgo.h.

Referenced by setPFEleParameters().

std::string PFAlgo::mvaWeightFileEleID_

private

Variables for PFElectrons.

Definition at line 360 of file PFAlgo.h.

Referenced by setPFEleParameters().

double PFAlgo::nSigmaECAL_

private

number of sigma to judge energy excess in ECAL

Definition at line 346 of file PFAlgo.h.

Referenced by operator<<(), processBlock(), and setParameters().

double PFAlgo::nSigmaHCAL_

private

number of sigma to judge energy excess in HCAL

Definition at line 349 of file PFAlgo.h.

Referenced by nSigmaHCAL(), operator<<(), and setParameters().

double PFAlgo::nSigmaTRACK_

private

Definition at line 404 of file PFAlgo.h.

Referenced by processBlock(), and setPFMuonAndFakeParameters().

unsigned int PFAlgo::nTrackIsoForEgammaSC_

private

Definition at line 375 of file PFAlgo.h.

Referenced by setPFEleParameters().

std::auto_ptr< reco::PFCandidateCollection > PFAlgo::pfAddedMuonCandidates_

protected

the collection of added muon candidates

Definition at line 311 of file PFAlgo.h.

Referenced by postMuonCleaning(), reconstructParticles(), and transferAddedMuonCandidates().

std::auto_ptr< reco::PFCandidateCollection > PFAlgo::pfCandidates_

protected

Definition at line 293 of file PFAlgo.h.

Referenced by checkCleaning(), operator<<(), pfCandidates(), postCleaning(), postMuonCleaning(), processBlock(), reconstructCluster(), reconstructParticles(), reconstructTrack(), setElectronExtraRef(), setPhotonExtraRef(), and transferCandidates().

std::auto_ptr< reco::PFCandidateCollection > PFAlgo::pfCleanedCandidates_

protected

Definition at line 299 of file PFAlgo.h.

Referenced by postCleaning(), reconstructParticles(), and transferCleanedCandidates().

std::auto_ptr< reco::PFCandidateCollection > PFAlgo::pfCleanedTrackerAndGlobalMuonCandidates_

protected

the collection of tracker/global cleaned muon candidates

Definition at line 303 of file PFAlgo.h.

Referenced by postMuonCleaning(), reconstructParticles(), and transferCleanedTrackerAndGlobalMuonCandidates().

std::auto_ptr< reco::PFCandidateCollection > PFAlgo::pfCosmicsMuonCleanedCandidates_

protected

the collection of cosmics cleaned muon candidates

Definition at line 301 of file PFAlgo.h.

Referenced by postMuonCleaning(), reconstructParticles(), and transferCosmicsMuonCleanedCandidates().

PFElectronAlgo* PFAlgo::pfele_

private

Definition at line 376 of file PFAlgo.h.

Referenced by processBlock(), setEGElectronCollection(), setPFEleParameters(), and ~PFAlgo().

std::auto_ptr< reco::PFCandidateCollection > PFAlgo::pfElectronCandidates_

protected

the unfiltered electron collection

Definition at line 295 of file PFAlgo.h.

Referenced by processBlock(), reconstructParticles(), setElectronExtraRef(), and transferElectronCandidates().

reco::PFCandidateElectronExtraCollection PFAlgo::pfElectronExtra_

protected

the unfiltered electron collection

Definition at line 314 of file PFAlgo.h.

Referenced by processBlock(), reconstructParticles(), setElectronExtraRef(), and transferElectronExtra().

std::auto_ptr< reco::PFCandidateCollection > PFAlgo::pfFakeMuonCleanedCandidates_

protected

the collection of fake cleaned muon candidates

Definition at line 305 of file PFAlgo.h.

Referenced by postMuonCleaning(), reconstructParticles(), and transferFakeMuonCleanedCandidates().

PFPhotonAlgo* PFAlgo::pfpho_

private

Definition at line 377 of file PFAlgo.h.

Referenced by processBlock(), setPFPhotonParameters(), setPFPhotonRegWeights(), setPFVertexParameters(), and ~PFAlgo().

std::auto_ptr< reco::PFCandidateCollection > PFAlgo::pfPhotonCandidates_

protected

the unfiltered photon collection

Definition at line 297 of file PFAlgo.h.

Referenced by processBlock(), and reconstructParticles().

reco::PFCandidatePhotonExtraCollection PFAlgo::pfPhotonExtra_

protected

the extra photon collection

Definition at line 316 of file PFAlgo.h.

Referenced by processBlock(), reconstructParticles(), setPhotonExtraRef(), and transferPhotonExtra().

std::auto_ptr< reco::PFCandidateCollection > PFAlgo::pfPunchThroughHadronCleanedCandidates_

protected

the collection of punch-through cleaned neutral hadron candidates

Definition at line 309 of file PFAlgo.h.

Referenced by postMuonCleaning(), reconstructParticles(), and transferPunchThroughHadronCleanedCandidates().

std::auto_ptr< reco::PFCandidateCollection > PFAlgo::pfPunchThroughMuonCleanedCandidates_

protected

the collection of punch-through cleaned muon candidates

Definition at line 307 of file PFAlgo.h.

Referenced by postMuonCleaning(), reconstructParticles(), and transferPunchThroughMuonCleanedCandidates().

bool PFAlgo::postHFCleaning_

private

Definition at line 409 of file PFAlgo.h.

Referenced by postCleaning(), and setPostHFCleaningParameters().

bool PFAlgo::postMuonCleaning_

private

Definition at line 410 of file PFAlgo.h.

reco::Vertex PFAlgo::primaryVertex_

private

Definition at line 420 of file PFAlgo.h.

Referenced by postMuonCleaning(), processBlock(), reconstructCluster(), setPFPhotonParameters(), and setPFVertexParameters().

double PFAlgo::ptError_

private

Definition at line 405 of file PFAlgo.h.

Referenced by processBlock(), and setPFMuonAndFakeParameters().

bool PFAlgo::rejectTracks_Bad_

private

Flags to use the protection against fakes and not reconstructed displaced vertices

Definition at line 384 of file PFAlgo.h.

Referenced by processBlock(), and setDisplacedVerticesParameters().

bool PFAlgo::rejectTracks_Step45_

private

Definition at line 385 of file PFAlgo.h.

Referenced by processBlock(), and setDisplacedVerticesParameters().

std::vector<double> PFAlgo::setchi2Values_

private

Definition at line 361 of file PFAlgo.h.

double PFAlgo::sumEtEcalIsoForEgammaSC_barrel_

private

Definition at line 369 of file PFAlgo.h.

Referenced by setPFEleParameters().

double PFAlgo::sumEtEcalIsoForEgammaSC_endcap_

private

Definition at line 370 of file PFAlgo.h.

Referenced by setPFEleParameters().

double PFAlgo::sumPtTrackIsoForEgammaSC_barrel_

private

Definition at line 372 of file PFAlgo.h.

Referenced by setPFEleParameters().

double PFAlgo::sumPtTrackIsoForEgammaSC_endcap_

private

Definition at line 373 of file PFAlgo.h.

Referenced by setPFEleParameters().

boost::shared_ptr<PFEnergyCalibrationHF> PFAlgo::thepfEnergyCalibrationHF_

private

Definition at line 352 of file PFAlgo.h.

Referenced by processBlock(), and setParameters().

boost::shared_ptr<PFSCEnergyCalibration> PFAlgo::thePFSCEnergyCalibration_

private

Definition at line 353 of file PFAlgo.h.

Referenced by setPFEleParameters().

double PFAlgo::useBestMuonTrack_

private

Definition at line 417 of file PFAlgo.h.

Referenced by postMuonCleaning(), processBlock(), reconstructTrack(), and setPFMuonAndFakeParameters().

bool PFAlgo::useEGammaSupercluster_

private

Definition at line 368 of file PFAlgo.h.

Referenced by setPFEleParameters().

bool PFAlgo::useEGElectrons_

private

Definition at line 367 of file PFAlgo.h.

Referenced by setEGElectronCollection(), and setPFEleParameters().

bool PFAlgo::useHO_

private

Definition at line 355 of file PFAlgo.h.

Referenced by processBlock(), and setHOTag().

bool PFAlgo::usePFConversions_

private

Definition at line 388 of file PFAlgo.h.

Referenced by processBlock(), and setDisplacedVerticesParameters().

bool PFAlgo::usePFDecays_

private

Definition at line 389 of file PFAlgo.h.

Referenced by isFromSecInt(), and setDisplacedVerticesParameters().

bool PFAlgo::usePFElectrons_

private

Definition at line 363 of file PFAlgo.h.

Referenced by processBlock(), setElectronExtraRef(), setPFEleParameters(), and ~PFAlgo().

bool PFAlgo::usePFMuonMomAssign_

private

Definition at line 380 of file PFAlgo.h.

Referenced by reconstructTrack(), and setPFMuonAndFakeParameters().

bool PFAlgo::usePFNuclearInteractions_

private

Definition at line 387 of file PFAlgo.h.

Referenced by isFromSecInt(), and setDisplacedVerticesParameters().

bool PFAlgo::usePFPhotons_

private

Definition at line 364 of file PFAlgo.h.

Referenced by processBlock(), setPFPhotonParameters(), setPFVertexParameters(), setPhotonExtraRef(), and ~PFAlgo().

bool PFAlgo::usePFSCEleCalib_

private

Definition at line 366 of file PFAlgo.h.

Referenced by setPFEleParameters().

bool PFAlgo::useVertices_

private

Definition at line 421 of file PFAlgo.h.

Referenced by reconstructCluster(), setPFPhotonParameters(), and setPFVertexParameters().

Public Member Functions

Protected Member Functions

Protected Attributes

Private Member Functions

Private Attributes

Friends

Detailed Description

Constructor & Destructor Documentation

Member Function Documentation

Friends And Related Function Documentation

Member Data Documentation