|
|
|
Particle Flow Algorithm test bench for the electron team. More...
#include <PFAlgoTestBenchElectrons.h>
Public Member Functions | |
| PFAlgoTestBenchElectrons () | |
| constructor | |
| virtual | ~PFAlgoTestBenchElectrons () |
| destructor | |
Protected Member Functions | |
| virtual void | processBlock (const reco::PFBlockRef &blockref, std::list< reco::PFBlockRef > &hcalBlockRefs, std::list< reco::PFBlockRef > &ecalBlockRefs) |
Particle Flow Algorithm test bench for the electron team.
Definition at line 18 of file PFAlgoTestBenchElectrons.h.
| PFAlgoTestBenchElectrons::PFAlgoTestBenchElectrons | ( | ) | [inline] |
| virtual PFAlgoTestBenchElectrons::~PFAlgoTestBenchElectrons | ( | ) | [inline, virtual] |
| void PFAlgoTestBenchElectrons::processBlock | ( | const reco::PFBlockRef & | blockref, |
| std::list< reco::PFBlockRef > & | hcalBlockRefs, | ||
| std::list< reco::PFBlockRef > & | ecalBlockRefs | ||
| ) | [protected, virtual] |
process one block. can be reimplemented in more sophisticated algorithms
Reimplemented from PFAlgo.
Definition at line 451 of file PFAlgo.cc.
References a, abs, PFAlgo::associatePSClusters(), b, edm::OwnVector< T, P >::begin(), Association::block, alignmentValidation::c1, PFAlgo::calibration_, dtNoiseDBValidation_cfg::cerr, CastorDataFrameFilter_impl::check(), gather_cfg::cout, PFAlgo::debug_, PFAlgo::dptRel_DispVtx_, ECAL, reco::PFBlockElement::ECAL, asciidump::elements, reco::PFCandidate::elementsInBlocks(), reco::LeafCandidate::energy(), eta(), PFAlgo::factors45_, first, PFElectronAlgo::getAllElectronCandidates(), PFElectronAlgo::getElectronCandidates(), PFElectronAlgo::getElectronExtra(), reco::PFBlockElement::GSF, HCAL, reco::PFBlockElement::HCAL, PFLayer::HF_EM, PFLayer::HF_HAD, reco::TrackBase::highPurity, reco::PFBlockElement::HO, i, getHLTprescales::index, PFElectronAlgo::isElectronValidCandidate(), PFAlgo::isFromSecInt(), PFMuonAlgo::isGlobalLooseMuon(), PFMuonAlgo::isIsolatedMuon(), PFMuonAlgo::isLooseMuon(), reco::isMuon(), edm::Ref< C, T, F >::isNull(), PFPhotonAlgo::isPhotonValidCandidate(), j, reco::PFBlock::LINKTEST_ALL, max(), min, RPCpg::mu, PFAlgo::muonECAL_, PFAlgo::muonHCAL_, PFAlgo::muonHO_, reco::btau::neutralEnergy, PFAlgo::neutralHadronEnergyResolution(), PFAlgo::nSigmaECAL_, PFAlgo::nSigmaHCAL(), PFAlgo::nSigmaTRACK_, convertSQLiteXML::ok, AlCaHLTBitMon_ParallelJobs::p, PFAlgo::pfCandidates_, PFAlgo::pfele_, PFAlgo::pfElectronCandidates_, PFAlgo::pfElectronExtra_, PFAlgo::pfpho_, PFAlgo::pfPhotonCandidates_, PFAlgo::pfPhotonExtra_, PFAlgo::primaryVertex_, PFMuonAlgo::printMuonProperties(), reco::PFBlockElement::PS1, reco::PFBlockElement::PS2, PFAlgo::ptError_, edm::OwnVector< T, P >::push_back(), PFAlgo::reconstructCluster(), PFAlgo::reconstructTrack(), PFAlgo::rejectTracks_Bad_, PFAlgo::rejectTracks_Step45_, edm::second(), reco::PFCandidate::setEcalEnergy(), reco::PFCandidate::setHcalEnergy(), reco::PFCandidate::setHoEnergy(), reco::PFBlockElementTrack::setPositionAtECALEntrance(), edm::OwnVector< T, P >::size(), createPayload::skip, mathSSE::sqrt(), reco::PFBlockElement::T_FROM_GAMMACONV, PFAlgo::thepfEnergyCalibrationHF_, reco::PFBlockElement::TRACK, reco::btau::trackMomentum, reco::PFBlockElementTrack::trackType(), funct::true, PFAlgo::useBestMuonTrack_, PFAlgo::useHO_, PFAlgo::usePFConversions_, PFAlgo::usePFElectrons_, PFAlgo::usePFPhotons_, X, x, and Gflash::Z.
{
// 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
1.7.3