|
|
|
#include <PFPhotonAlgo.h>
Public Member Functions | |
| bool | isPhotonValidCandidate (const reco::PFBlockRef &blockRef, std::vector< bool > &active, std::auto_ptr< reco::PFCandidateCollection > &pfPhotonCandidates, std::auto_ptr< reco::PFCandidateCollection > &) |
| PFPhotonAlgo (std::string mvaweightfile, double mvaConvCut, const reco::Vertex &primary, const boost::shared_ptr< PFEnergyCalibration > &thePFEnergyCalibration) | |
| ~PFPhotonAlgo () | |
Private Types | |
| enum | verbosityLevel { Silent, Summary, Chatty } |
Private Member Functions | |
| bool | EvaluateSingleLegMVA (const reco::PFBlockRef &blockref, const reco::Vertex &primaryvtx, unsigned int track_index) |
| void | RunPFPhoton (const reco::PFBlockRef &blockRef, std::vector< bool > &active, std::auto_ptr< reco::PFCandidateCollection > &pfPhotonCandidates) |
Private Attributes | |
| float | chi2 |
| float | del_phi |
| float | EoverPt |
| float | HoverPt |
| bool | isvalid_ |
| double | MVACUT |
| double | mvaValue |
| float | nlayers |
| float | nlost |
| reco::Vertex | primaryVertex_ |
| float | STIP |
| boost::shared_ptr < PFEnergyCalibration > | thePFEnergyCalibration_ |
| TMVA::Reader * | tmvaReader_ |
| float | track_pt |
| verbosityLevel | verbosityLevel_ |
Definition at line 24 of file PFPhotonAlgo.h.
enum PFPhotonAlgo::verbosityLevel [private] |
Definition at line 56 of file PFPhotonAlgo.h.
| PFPhotonAlgo::PFPhotonAlgo | ( | std::string | mvaweightfile, |
| double | mvaConvCut, | ||
| const reco::Vertex & | primary, | ||
| const boost::shared_ptr< PFEnergyCalibration > & | thePFEnergyCalibration | ||
| ) |
Definition at line 23 of file PFPhotonAlgo.cc.
References chi2, del_phi, EoverPt, HoverPt, nlayers, nlost, primaryVertex_, STIP, tmvaReader_, and track_pt.
: isvalid_(false), verbosityLevel_(Silent), MVACUT(mvaConvCut), thePFEnergyCalibration_(thePFEnergyCalibration) { primaryVertex_=primary; //Book MVA tmvaReader_ = new TMVA::Reader("!Color:Silent"); tmvaReader_->AddVariable("del_phi",&del_phi); tmvaReader_->AddVariable("nlayers", &nlayers); tmvaReader_->AddVariable("chi2",&chi2); tmvaReader_->AddVariable("EoverPt",&EoverPt); tmvaReader_->AddVariable("HoverPt",&HoverPt); tmvaReader_->AddVariable("track_pt", &track_pt); tmvaReader_->AddVariable("STIP",&STIP); tmvaReader_->AddVariable("nlost", &nlost); tmvaReader_->BookMVA("BDT",mvaweightfile.c_str()); }
| PFPhotonAlgo::~PFPhotonAlgo | ( | ) | [inline] |
| bool PFPhotonAlgo::EvaluateSingleLegMVA | ( | const reco::PFBlockRef & | blockref, |
| const reco::Vertex & | primaryvtx, | ||
| unsigned int | track_index | ||
| ) | [private] |
Definition at line 620 of file PFPhotonAlgo.cc.
References reco::PFBlock::associatedElements(), createPayload::block, chi2, del_phi, Geom::deltaPhi(), reco::PFBlockElement::ECAL, reco::PFBlock::elements(), asciidump::elements, EoverPt, reco::PFBlockElement::HCAL, HoverPt, reco::PFBlock::linkData(), reco::PFBlock::LINKTEST_ALL, MVACUT, mvaValue, nlayers, nlost, PV3DBase< T, PVType, FrameType >::phi(), colinearityKinematic::Phi, STIP, tmvaReader_, track_pt, reco::Vertex::x(), X, reco::Vertex::y(), and reco::Vertex::z().
Referenced by RunPFPhoton().
{
bool convtkfound=false;
const reco::PFBlock& block = *blockref;
const edm::OwnVector< reco::PFBlockElement >& elements = block.elements();
//use this to store linkdata in the associatedElements function below
PFBlock::LinkData linkData = block.linkData();
//calculate MVA Variables
chi2=elements[track_index].trackRef()->chi2()/elements[track_index].trackRef()->ndof();
nlost=elements[track_index].trackRef()->trackerExpectedHitsInner().numberOfLostHits();
nlayers=elements[track_index].trackRef()->hitPattern().trackerLayersWithMeasurement();
track_pt=elements[track_index].trackRef()->pt();
STIP=elements[track_index].trackRefPF()->STIP();
float linked_e=0;
float linked_h=0;
std::multimap<double, unsigned int> ecalAssoTrack;
block.associatedElements( track_index,linkData,
ecalAssoTrack,
reco::PFBlockElement::ECAL,
reco::PFBlock::LINKTEST_ALL );
std::multimap<double, unsigned int> hcalAssoTrack;
block.associatedElements( track_index,linkData,
hcalAssoTrack,
reco::PFBlockElement::HCAL,
reco::PFBlock::LINKTEST_ALL );
if(ecalAssoTrack.size() > 0) {
for(std::multimap<double, unsigned int>::iterator itecal = ecalAssoTrack.begin();
itecal != ecalAssoTrack.end(); ++itecal) {
linked_e=linked_e+elements[itecal->second].clusterRef()->energy();
}
}
if(hcalAssoTrack.size() > 0) {
for(std::multimap<double, unsigned int>::iterator ithcal = hcalAssoTrack.begin();
ithcal != hcalAssoTrack.end(); ++ithcal) {
linked_h=linked_h+elements[ithcal->second].clusterRef()->energy();
}
}
EoverPt=linked_e/elements[track_index].trackRef()->pt();
HoverPt=linked_h/elements[track_index].trackRef()->pt();
GlobalVector rvtx(elements[track_index].trackRef()->innerPosition().X()-primaryvtx.x(),
elements[track_index].trackRef()->innerPosition().Y()-primaryvtx.y(),
elements[track_index].trackRef()->innerPosition().Z()-primaryvtx.z());
double vtx_phi=rvtx.phi();
//delta Phi between conversion vertex and track
del_phi=fabs(deltaPhi(vtx_phi, elements[track_index].trackRef()->innerMomentum().Phi()));
mvaValue = tmvaReader_->EvaluateMVA("BDT");
if(mvaValue > MVACUT)convtkfound=true;
return convtkfound;
}
| bool PFPhotonAlgo::isPhotonValidCandidate | ( | const reco::PFBlockRef & | blockRef, |
| std::vector< bool > & | active, | ||
| std::auto_ptr< reco::PFCandidateCollection > & | pfPhotonCandidates, | ||
| std::auto_ptr< reco::PFCandidateCollection > & | |||
| ) | [inline] |
Definition at line 37 of file PFPhotonAlgo.h.
References isvalid_, and RunPFPhoton().
Referenced by PFAlgoTestBenchElectrons::processBlock().
{
isvalid_=false;
// RunPFPhoton has to set isvalid_ to TRUE in case it finds a valid candidate
// ... TODO: maybe can be replaced by checking for the size of the CandCollection.....
RunPFPhoton(blockRef,
active,
pfPhotonCandidates);
return isvalid_;
};
| void PFPhotonAlgo::RunPFPhoton | ( | const reco::PFBlockRef & | blockRef, |
| std::vector< bool > & | active, | ||
| std::auto_ptr< reco::PFCandidateCollection > & | pfPhotonCandidates | ||
| ) | [private] |
Definition at line 46 of file PFPhotonAlgo.cc.
References reco::PFCandidate::addElementInBlock(), edm::OwnVector< T, P >::begin(), Chatty, convBrem_cff::convTracks, ECAL, asciidump::elements, edm::OwnVector< T, P >::end(), EvaluateSingleLegMVA(), reco::PFBlockElementSuperCluster::fromPhoton(), reco::PFCandidate::gamma, reco::PFCandidate::GAMMA_TO_GAMMACONV, reco::PFBlockElement::HCAL, i, isvalid_, reco::PFBlock::LINKTEST_ALL, mvaValue, primaryVertex_, reco::PFBlockElement::PS1, reco::PFBlockElement::PS2, edm::OwnVector< T, P >::push_back(), reco::PFBlockElement::SC, reco::PFCandidate::set_mva_nothing_gamma(), reco::PFCandidate::setEcalEnergy(), reco::PFCandidate::setFlag(), reco::PFCandidate::setHcalEnergy(), reco::PFCandidate::setPs1Energy(), reco::PFCandidate::setPs2Energy(), reco::PFCandidate::setSuperClusterRef(), edm::OwnVector< T, P >::size(), mathSSE::sqrt(), reco::PFBlockElementSuperCluster::superClusterRef(), reco::PFBlockElement::T_FROM_GAMMACONV, thePFEnergyCalibration_, reco::PFBlockElement::TRACK, reco::PFBlockElementTrack::trackType(), and verbosityLevel_.
Referenced by isPhotonValidCandidate().
{
//std::cout<<" calling RunPFPhoton "<<std::endl;
/* For now we construct the PhotonCandidate simply from
a) adding the CORRECTED energies of each participating ECAL cluster
b) build the energy-weighted direction for the Photon
*/
// define how much is printed out for debugging.
// ... will be setable via CFG file parameter
verbosityLevel_ = Chatty; // Chatty mode.
// loop over all elements in the Block
const edm::OwnVector< reco::PFBlockElement >& elements = blockRef->elements();
edm::OwnVector< reco::PFBlockElement >::const_iterator ele = elements.begin();
std::vector<bool>::const_iterator actIter = active.begin();
PFBlock::LinkData linkData = blockRef->linkData();
bool isActive = true;
if(elements.size() != active.size()) {
// throw excpetion...
//std::cout<<" WARNING: Size of collection and active-vectro don't agree!"<<std::endl;
return;
}
// local vecotr to keep track of the indices of the 'elements' for the Photon candidate
// once we decide to keep the candidate, the 'active' entriesd for them must be set to false
std::vector<unsigned int> elemsToLock;
elemsToLock.resize(0);
for( ; ele != elements.end(); ++ele, ++actIter ) {
// if it's not a SuperCluster, go to the next element
if( !( ele->type() == reco::PFBlockElement::SC ) ) continue;
// Photon kienmatics, will be updated for each identified participating element
float photonEnergy_ = 0.;
float photonX_ = 0.;
float photonY_ = 0.;
float photonZ_ = 0.;
float RawEcalEne = 0.;
// Total pre-shower energy
float ps1TotEne = 0.;
float ps2TotEne = 0.;
bool hasConvTrack=false;
bool hasSingleleg=false;
std::vector<unsigned int> AddClusters(0);
std::vector<unsigned int> IsoTracks(0);
std::multimap<unsigned int, unsigned int>ClusterAddPS1;
std::multimap<unsigned int, unsigned int>ClusterAddPS2;
isActive = *(actIter);
//cout << " Found a SuperCluster. Energy " ;
const reco::PFBlockElementSuperCluster *sc = dynamic_cast<const reco::PFBlockElementSuperCluster*>(&(*ele));
//std::cout << sc->superClusterRef()->energy () << " Track/Ecal/Hcal Iso " << sc->trackIso()<< " " << sc->ecalIso() ;
//std::cout << " " << sc->hcalIso() <<std::endl;
if (!(sc->fromPhoton()))continue;
// check the status of the SC Element...
// ..... I understand it should *always* be active, since PFElectronAlgo does not touch this (yet?) RISHI: YES
if( !isActive ) {
//std::cout<<" SuperCluster is NOT active.... "<<std::endl;
continue;
}
elemsToLock.push_back(ele-elements.begin()); //add SC to elements to lock
// loop over its constituent ECAL cluster
std::multimap<double, unsigned int> ecalAssoPFClusters;
blockRef->associatedElements( ele-elements.begin(),
linkData,
ecalAssoPFClusters,
reco::PFBlockElement::ECAL,
reco::PFBlock::LINKTEST_ALL );
// loop over the ECAL clusters linked to the iEle
if( ! ecalAssoPFClusters.size() ) {
// This SC element has NO ECAL elements asigned... *SHOULD NOT HAPPEN*
//std::cout<<" Found SC element with no ECAL assigned "<<std::endl;
continue;
}
// This is basically CASE 2
// .... we loop over all ECAL cluster linked to each other by this SC
for(std::multimap<double, unsigned int>::iterator itecal = ecalAssoPFClusters.begin();
itecal != ecalAssoPFClusters.end(); ++itecal) {
// to get the reference to the PF clusters, this is needed.
reco::PFClusterRef clusterRef = elements[itecal->second].clusterRef();
// from the clusterRef get the energy, direction, etc
// float ClustRawEnergy = clusterRef->energy();
// float ClustEta = clusterRef->position().eta();
// float ClustPhi = clusterRef->position().phi();
// initialize the vectors for the PS energies
vector<double> ps1Ene(0);
vector<double> ps2Ene(0);
double ps1=0;
double ps2=0;
hasSingleleg=false;
hasConvTrack=false;
/*
cout << " My cluster index " << itecal->second
<< " energy " << ClustRawEnergy
<< " eta " << ClustEta
<< " phi " << ClustPhi << endl;
*/
// check if this ECAL element is still active (could have been eaten by PFElectronAlgo)
// ......for now we give the PFElectron Algo *ALWAYS* Shot-Gun on the ECAL elements to the PFElectronAlgo
if( !( active[itecal->second] ) ) {
//std::cout<< " .... this ECAL element is NOT active anymore. Is skipped. "<<std::endl;
continue;
}
// ------------------------------------------------------------------------------------------
// TODO: do some tests on the ECAL cluster itself, deciding to use it or not for the Photons
// ..... ??? Do we need this?
if ( false ) {
// Check if there are a large number tracks that do not pass pre-ID around this ECAL cluster
bool useIt = true;
int mva_reject=0;
bool isClosest=false;
std::multimap<double, unsigned int> Trackscheck;
blockRef->associatedElements( itecal->second,
linkData,
Trackscheck,
reco::PFBlockElement::TRACK,
reco::PFBlock::LINKTEST_ALL);
for(std::multimap<double, unsigned int>::iterator track = Trackscheck.begin();
track != Trackscheck.end(); ++track) {
// first check if is it's still active
if( ! (active[track->second]) ) continue;
hasSingleleg=EvaluateSingleLegMVA(blockRef, primaryVertex_, track->second);
//check if it is the closest linked track
std::multimap<double, unsigned int> closecheck;
blockRef->associatedElements(track->second,
linkData,
closecheck,
reco::PFBlockElement::ECAL,
reco::PFBlock::LINKTEST_ALL);
if(closecheck.begin()->second ==itecal->second)isClosest=true;
if(!hasSingleleg)mva_reject++;
}
if(mva_reject>0 && isClosest)useIt=false;
//if(mva_reject==1 && isClosest)useIt=false;
if( !useIt ) continue; // Go to next ECAL cluster within SC
}
// ------------------------------------------------------------------------------------------
// We decided to keep the ECAL cluster for this Photon Candidate ...
elemsToLock.push_back(itecal->second);
// look for PS in this Block linked to this ECAL cluster
std::multimap<double, unsigned int> PS1Elems;
std::multimap<double, unsigned int> PS2Elems;
//PS Layer 1 linked to ECAL cluster
blockRef->associatedElements( itecal->second,
linkData,
PS1Elems,
reco::PFBlockElement::PS1,
reco::PFBlock::LINKTEST_ALL );
//PS Layer 2 linked to the ECAL cluster
blockRef->associatedElements( itecal->second,
linkData,
PS2Elems,
reco::PFBlockElement::PS2,
reco::PFBlock::LINKTEST_ALL );
// loop over all PS1 and compute energy
for(std::multimap<double, unsigned int>::iterator iteps = PS1Elems.begin();
iteps != PS1Elems.end(); ++iteps) {
// first chekc if it's still active
if( !(active[iteps->second]) ) continue;
//Check if this PS1 is not closer to another ECAL cluster in this Block
std::multimap<double, unsigned int> ECALPS1check;
blockRef->associatedElements( iteps->second,
linkData,
ECALPS1check,
reco::PFBlockElement::ECAL,
reco::PFBlock::LINKTEST_ALL );
if(itecal->second==ECALPS1check.begin()->second)//then it is closest linked
{
reco::PFClusterRef ps1ClusterRef = elements[iteps->second].clusterRef();
ps1Ene.push_back( ps1ClusterRef->energy() );
ps1=ps1+ps1ClusterRef->energy(); //add to total PS1
// incativate this PS1 Element
elemsToLock.push_back(iteps->second);
}
}
for(std::multimap<double, unsigned int>::iterator iteps = PS2Elems.begin();
iteps != PS2Elems.end(); ++iteps) {
// first chekc if it's still active
if( !(active[iteps->second]) ) continue;
// Check if this PS2 is not closer to another ECAL cluster in this Block:
std::multimap<double, unsigned int> ECALPS2check;
blockRef->associatedElements( iteps->second,
linkData,
ECALPS2check,
reco::PFBlockElement::ECAL,
reco::PFBlock::LINKTEST_ALL );
if(itecal->second==ECALPS2check.begin()->second)//is closest linked
{
reco::PFClusterRef ps2ClusterRef = elements[iteps->second].clusterRef();
ps2Ene.push_back( ps2ClusterRef->energy() );
ps2=ps2ClusterRef->energy()+ps2; //add to total PS2
// incativate this PS2 Element
elemsToLock.push_back(iteps->second);
}
}
// loop over the HCAL Clusters linked to the ECAL cluster (CASE 6)
std::multimap<double, unsigned int> hcalElems;
blockRef->associatedElements( itecal->second,linkData,
hcalElems,
reco::PFBlockElement::HCAL,
reco::PFBlock::LINKTEST_ALL );
for(std::multimap<double, unsigned int>::iterator ithcal = hcalElems.begin();
ithcal != hcalElems.end(); ++ithcal) {
if ( ! (active[ithcal->second] ) ) continue; // HCAL Cluster already used....
// TODO: Decide if this HCAL cluster is to be used
// .... based on some Physics
// .... To we need to check if it's closer to any other ECAL/TRACK?
bool useHcal = false;
if ( !useHcal ) continue;
//not locked
//elemsToLock.push_back(ithcal->second);
}
// This is entry point for CASE 3.
// .... we loop over all Tracks linked to this ECAL and check if it's labeled as conversion
// This is the part for looping over all 'Conversion' Tracks
std::multimap<double, unsigned int> convTracks;
blockRef->associatedElements( itecal->second,
linkData,
convTracks,
reco::PFBlockElement::TRACK,
reco::PFBlock::LINKTEST_ALL);
for(std::multimap<double, unsigned int>::iterator track = convTracks.begin();
track != convTracks.end(); ++track) {
// first check if is it's still active
if( ! (active[track->second]) ) continue;
// check if it's a CONV track
const reco::PFBlockElementTrack * trackRef = dynamic_cast<const reco::PFBlockElementTrack*>((&elements[track->second]));
//Check if track is a Single leg from a Conversion
mvaValue=-999;
hasSingleleg=EvaluateSingleLegMVA(blockRef, primaryVertex_, track->second);
//If it is not then it will be used to check Track Isolation at the end
if(!hasSingleleg)
{
bool included=false;
//check if this track is already included in the vector so it is linked to an ECAL cluster that is already examined
for(unsigned int i=0; i<IsoTracks.size(); i++)
{if(IsoTracks[i]==track->second)included=true;}
if(!included)IsoTracks.push_back(track->second);
}
//For now only Pre-ID tracks that are not already identified as Conversions
if(hasSingleleg &&!(trackRef->trackType(reco::PFBlockElement::T_FROM_GAMMACONV)))
{
elemsToLock.push_back(track->second);
//find all the clusters linked to this track
std::multimap<double, unsigned int> moreClusters;
blockRef->associatedElements( track->second,
linkData,
moreClusters,
reco::PFBlockElement::ECAL,
reco::PFBlock::LINKTEST_ALL);
float p_in=sqrt(elements[track->second].trackRef()->innerMomentum().x() * elements[track->second].trackRef()->innerMomentum().x() +
elements[track->second].trackRef()->innerMomentum().y()*elements[track->second].trackRef()->innerMomentum().y()+
elements[track->second].trackRef()->innerMomentum().z()*elements[track->second].trackRef()->innerMomentum().z());
float linked_E=0;
for(std::multimap<double, unsigned int>::iterator clust = moreClusters.begin();
clust != moreClusters.end(); ++clust)
{
if(!active[clust->second])continue;
//running sum of linked energy
linked_E=linked_E+elements[clust->second].clusterRef()->energy();
//prevent too much energy from being added
if(linked_E/p_in>1.5)break;
bool included=false;
//check if these ecal clusters are already included with the supercluster
for(std::multimap<double, unsigned int>::iterator cluscheck = ecalAssoPFClusters.begin();
cluscheck != ecalAssoPFClusters.end(); ++cluscheck)
{
if(cluscheck->second==clust->second)included=true;
}
if(!included)AddClusters.push_back(clust->second);//Add to a container of clusters to be Added to the Photon candidate
}
}
// Possibly need to be more smart about them (CASE 5)
// .... for now we simply skip non id'ed tracks
if( ! (trackRef->trackType(reco::PFBlockElement::T_FROM_GAMMACONV) ) ) continue;
hasConvTrack=true;
elemsToLock.push_back(track->second);
//again look at the clusters linked to this track
std::multimap<double, unsigned int> moreClusters;
blockRef->associatedElements( track->second,
linkData,
moreClusters,
reco::PFBlockElement::ECAL,
reco::PFBlock::LINKTEST_ALL);
float p_in=sqrt(elements[track->second].trackRef()->innerMomentum().x() * elements[track->second].trackRef()->innerMomentum().x() +
elements[track->second].trackRef()->innerMomentum().y()*elements[track->second].trackRef()->innerMomentum().y()+
elements[track->second].trackRef()->innerMomentum().z()*elements[track->second].trackRef()->innerMomentum().z());
float linked_E=0;
for(std::multimap<double, unsigned int>::iterator clust = moreClusters.begin();
clust != moreClusters.end(); ++clust)
{
if(!active[clust->second])continue;
linked_E=linked_E+elements[clust->second].clusterRef()->energy();
if(linked_E/p_in>1.5)break;
bool included=false;
for(std::multimap<double, unsigned int>::iterator cluscheck = ecalAssoPFClusters.begin();
cluscheck != ecalAssoPFClusters.end(); ++cluscheck)
{
if(cluscheck->second==clust->second)included=true;
}
if(!included)AddClusters.push_back(clust->second);//again only add if it is not already included with the supercluster
}
// we need to check for other TRACKS linked to this conversion track, that point possibly no an ECAL cluster not included in the SC
// .... This is basically CASE 4.
std::multimap<double, unsigned int> moreTracks;
blockRef->associatedElements( track->second,
linkData,
moreTracks,
reco::PFBlockElement::TRACK,
reco::PFBlock::LINKTEST_ALL);
for(std::multimap<double, unsigned int>::iterator track2 = moreTracks.begin();
track2 != moreTracks.end(); ++track2) {
// first check if is it's still active
if( ! (active[track2->second]) ) continue;
//skip over the 1st leg already found above
if(track->second==track2->second)continue;
// check if it's a CONV track
const reco::PFBlockElementTrack * track2Ref = dynamic_cast<const reco::PFBlockElementTrack*>((&elements[track2->second]));
if( ! (track2Ref->trackType(reco::PFBlockElement::T_FROM_GAMMACONV) ) ) continue; // Possibly need to be more smart about them (CASE 5)
elemsToLock.push_back(track2->second);
// so it's another active conversion track, that is in the Block and linked to the conversion track we already found
// find the ECAL cluster linked to it...
std::multimap<double, unsigned int> convEcal;
blockRef->associatedElements( track2->second,
linkData,
convEcal,
reco::PFBlockElement::ECAL,
reco::PFBlock::LINKTEST_ALL);
float p_in=sqrt(elements[track->second].trackRef()->innerMomentum().x()*elements[track->second].trackRef()->innerMomentum().x()+
elements[track->second].trackRef()->innerMomentum().y()*elements[track->second].trackRef()->innerMomentum().y()+
elements[track->second].trackRef()->innerMomentum().z()*elements[track->second].trackRef()->innerMomentum().z());
float linked_E=0;
for(std::multimap<double, unsigned int>::iterator itConvEcal = convEcal.begin();
itConvEcal != convEcal.end(); ++itConvEcal) {
if( ! (active[itConvEcal->second]) ) continue;
bool included=false;
for(std::multimap<double, unsigned int>::iterator cluscheck = ecalAssoPFClusters.begin();
cluscheck != ecalAssoPFClusters.end(); ++cluscheck)
{
if(cluscheck->second==itConvEcal->second)included=true;
}
linked_E=linked_E+elements[itConvEcal->second].clusterRef()->energy();
if(linked_E/p_in>1.5)break;
if(!included){AddClusters.push_back(itConvEcal->second);
}
// it's still active, so we have to add it.
// CAUTION: we don't care here if it's part of the SC or not, we include it anyways
// loop over the HCAL Clusters linked to the ECAL cluster (CASE 6)
std::multimap<double, unsigned int> hcalElems_conv;
blockRef->associatedElements( itecal->second,linkData,
hcalElems_conv,
reco::PFBlockElement::HCAL,
reco::PFBlock::LINKTEST_ALL );
for(std::multimap<double, unsigned int>::iterator ithcal2 = hcalElems_conv.begin();
ithcal2 != hcalElems_conv.end(); ++ithcal2) {
if ( ! (active[ithcal2->second] ) ) continue; // HCAL Cluster already used....
// TODO: Decide if this HCAL cluster is to be used
// .... based on some Physics
// .... To we need to check if it's closer to any other ECAL/TRACK?
bool useHcal = true;
if ( !useHcal ) continue;
//elemsToLock.push_back(ithcal2->second);
} // end of loop over HCAL clusters linked to the ECAL cluster from second CONVERSION leg
} // end of loop over ECALs linked to second T_FROM_GAMMACONV
} // end of loop over SECOND conversion leg
// TODO: Do we need to check separatly if there are HCAL cluster linked to the track?
} // end of loop over tracks
// Calibrate the Added ECAL energy
float addedCalibEne=0;
float addedRawEne=0;
std::vector<double>AddedPS1(0);
std::vector<double>AddedPS2(0);
double addedps1=0;
double addedps2=0;
for(unsigned int i=0; i<AddClusters.size(); i++)
{
std::multimap<double, unsigned int> PS1Elems_conv;
std::multimap<double, unsigned int> PS2Elems_conv;
blockRef->associatedElements(AddClusters[i],
linkData,
PS1Elems_conv,
reco::PFBlockElement::PS1,
reco::PFBlock::LINKTEST_ALL );
blockRef->associatedElements( AddClusters[i],
linkData,
PS2Elems_conv,
reco::PFBlockElement::PS2,
reco::PFBlock::LINKTEST_ALL );
for(std::multimap<double, unsigned int>::iterator iteps = PS1Elems_conv.begin();
iteps != PS1Elems_conv.end(); ++iteps)
{
if(!active[iteps->second])continue;
std::multimap<double, unsigned int> PS1Elems_check;
blockRef->associatedElements(iteps->second,
linkData,
PS1Elems_check,
reco::PFBlockElement::ECAL,
reco::PFBlock::LINKTEST_ALL );
if(PS1Elems_check.begin()->second==AddClusters[i])
{
reco::PFClusterRef ps1ClusterRef = elements[iteps->second].clusterRef();
AddedPS1.push_back(ps1ClusterRef->energy());
addedps1=addedps1+ps1ClusterRef->energy();
elemsToLock.push_back(iteps->second);
}
}
for(std::multimap<double, unsigned int>::iterator iteps = PS2Elems_conv.begin();
iteps != PS2Elems_conv.end(); ++iteps) {
if(!active[iteps->second])continue;
std::multimap<double, unsigned int> PS2Elems_check;
blockRef->associatedElements(iteps->second,
linkData,
PS2Elems_check,
reco::PFBlockElement::ECAL,
reco::PFBlock::LINKTEST_ALL );
if(PS2Elems_check.begin()->second==AddClusters[i])
{
reco::PFClusterRef ps2ClusterRef = elements[iteps->second].clusterRef();
AddedPS2.push_back(ps2ClusterRef->energy());
addedps2=addedps2+ps2ClusterRef->energy();
elemsToLock.push_back(iteps->second);
}
}
reco::PFClusterRef AddclusterRef = elements[AddClusters[i]].clusterRef();
addedRawEne = AddclusterRef->energy()+addedRawEne;
addedCalibEne = thePFEnergyCalibration_->energyEm(*AddclusterRef,AddedPS1,AddedPS2,false)+addedCalibEne;
AddedPS2.clear();
AddedPS1.clear();
elemsToLock.push_back(AddClusters[i]);
}
AddClusters.clear();
float EE = thePFEnergyCalibration_->energyEm(*clusterRef,ps1Ene,ps2Ene,false)+addedCalibEne;
//cout<<"Original Energy "<<EE<<"Added Energy "<<addedCalibEne<<endl;
photonEnergy_ += EE;
RawEcalEne += clusterRef->energy()+addedRawEne;
photonX_ += EE * clusterRef->position().X();
photonY_ += EE * clusterRef->position().Y();
photonZ_ += EE * clusterRef->position().Z();
ps1TotEne += ps1+addedps1;
ps2TotEne += ps2+addedps2;
} // end of loop over all ECAL cluster within this SC
// we've looped over all ECAL clusters, ready to generate PhotonCandidate
if( ! (photonEnergy_ > 0.) ) continue; // This SC is not a Photon Candidate
float sum_track_pt=0;
//Now check if there are tracks failing isolation outside of the Jurassic isolation region
for(unsigned int i=0; i<IsoTracks.size(); i++)sum_track_pt=sum_track_pt+elements[IsoTracks[i]].trackRef()->pt();
math::XYZVector photonPosition(photonX_,
photonY_,
photonZ_);
math::XYZVector photonDirection=photonPosition.Unit();
math::XYZTLorentzVector photonMomentum(photonEnergy_* photonDirection.X(),
photonEnergy_* photonDirection.Y(),
photonEnergy_* photonDirection.Z(),
photonEnergy_ );
if(sum_track_pt>(2+ 0.001* photonMomentum.pt()))
continue;//THIS SC is not a Photon it fails track Isolation
/*
std::cout<<" Created Photon with energy = "<<photonEnergy_<<std::endl;
std::cout<<" pT = "<<photonMomentum.pt()<<std::endl;
std::cout<<" RawEne = "<<RawEcalEne<<std::endl;
std::cout<<" E = "<<photonMomentum.e()<<std::endl;
std::cout<<" eta = "<<photonMomentum.eta()<<std::endl;
std::cout<<" TrackIsolation = "<< sum_track_pt <<std::endl;
*/
reco::PFCandidate photonCand(0,photonMomentum, reco::PFCandidate::gamma);
photonCand.setPs1Energy(ps1TotEne);
photonCand.setPs2Energy(ps2TotEne);
photonCand.setEcalEnergy(RawEcalEne,photonEnergy_);
photonCand.setHcalEnergy(0.,0.);
photonCand.set_mva_nothing_gamma(1.);
photonCand.setSuperClusterRef(sc->superClusterRef());
if(hasConvTrack || hasSingleleg)photonCand.setFlag( reco::PFCandidate::GAMMA_TO_GAMMACONV, true);
//photonCand.setPositionAtECALEntrance(math::XYZPointF(photonMom_.position()));
// set isvalid_ to TRUE since we've found at least one photon candidate
isvalid_ = true;
// push back the candidate into the collection ...
// ... and lock all elemts used
for(std::vector<unsigned int>::const_iterator it = elemsToLock.begin();
it != elemsToLock.end(); ++it)
{
if(active[*it])photonCand.addElementInBlock(blockRef,*it);
active[*it] = false;
}
pfCandidates->push_back(photonCand);
// ... and reset the vector
elemsToLock.resize(0);
hasConvTrack=false;
hasSingleleg=false;
} // end of loops over all elements in block
return;
}
float PFPhotonAlgo::chi2 [private] |
Definition at line 76 of file PFPhotonAlgo.h.
Referenced by EvaluateSingleLegMVA(), and PFPhotonAlgo().
float PFPhotonAlgo::del_phi [private] |
Definition at line 76 of file PFPhotonAlgo.h.
Referenced by EvaluateSingleLegMVA(), and PFPhotonAlgo().
float PFPhotonAlgo::EoverPt [private] |
Definition at line 76 of file PFPhotonAlgo.h.
Referenced by EvaluateSingleLegMVA(), and PFPhotonAlgo().
float PFPhotonAlgo::HoverPt [private] |
Definition at line 76 of file PFPhotonAlgo.h.
Referenced by EvaluateSingleLegMVA(), and PFPhotonAlgo().
bool PFPhotonAlgo::isvalid_ [private] |
Definition at line 63 of file PFPhotonAlgo.h.
Referenced by isPhotonValidCandidate(), and RunPFPhoton().
double PFPhotonAlgo::MVACUT [private] |
Definition at line 70 of file PFPhotonAlgo.h.
Referenced by EvaluateSingleLegMVA().
double PFPhotonAlgo::mvaValue [private] |
Definition at line 77 of file PFPhotonAlgo.h.
Referenced by EvaluateSingleLegMVA(), and RunPFPhoton().
float PFPhotonAlgo::nlayers [private] |
Definition at line 75 of file PFPhotonAlgo.h.
Referenced by EvaluateSingleLegMVA(), and PFPhotonAlgo().
float PFPhotonAlgo::nlost [private] |
Definition at line 75 of file PFPhotonAlgo.h.
Referenced by EvaluateSingleLegMVA(), and PFPhotonAlgo().
reco::Vertex PFPhotonAlgo::primaryVertex_ [private] |
Definition at line 71 of file PFPhotonAlgo.h.
Referenced by PFPhotonAlgo(), and RunPFPhoton().
float PFPhotonAlgo::STIP [private] |
Definition at line 76 of file PFPhotonAlgo.h.
Referenced by EvaluateSingleLegMVA(), and PFPhotonAlgo().
boost::shared_ptr<PFEnergyCalibration> PFPhotonAlgo::thePFEnergyCalibration_ [private] |
Definition at line 73 of file PFPhotonAlgo.h.
Referenced by RunPFPhoton().
TMVA::Reader* PFPhotonAlgo::tmvaReader_ [private] |
Definition at line 72 of file PFPhotonAlgo.h.
Referenced by EvaluateSingleLegMVA(), PFPhotonAlgo(), and ~PFPhotonAlgo().
float PFPhotonAlgo::track_pt [private] |
Definition at line 76 of file PFPhotonAlgo.h.
Referenced by EvaluateSingleLegMVA(), and PFPhotonAlgo().
verbosityLevel PFPhotonAlgo::verbosityLevel_ [private] |
Definition at line 64 of file PFPhotonAlgo.h.
Referenced by RunPFPhoton().
1.7.3