Skip to content

Commit

Permalink
Merge pull request #2490 from lgray/pfclusteralgo_stability
Browse files Browse the repository at this point in the history
Reco updates -- Protection to stabilize the behavior of PFClusterAlgo
  • Loading branch information
ktf committed Feb 22, 2014
2 parents 1910581 + f88a2d9 commit 4d93095
Show file tree
Hide file tree
Showing 5 changed files with 269 additions and 4 deletions.
Expand Up @@ -35,7 +35,7 @@
# n neighbours in PS
nNeighbours = cms.int32(8),
# sigma of the shower in preshower
showerSigma = cms.double(0.2),
showerSigma = cms.double(0.3),
# n crystals for position calculation in PS
posCalcNCrystal = cms.int32(-1),
# use cells with common corner to build topo-clusters
Expand Down
18 changes: 15 additions & 3 deletions RecoParticleFlow/PFClusterProducer/src/PFClusterAlgo.cc
Expand Up @@ -1122,6 +1122,13 @@ PFClusterAlgo::buildPFClusters( const std::vector< unsigned >& topocluster,

const reco::PFRecHit& rh = rechit( rhindex, rechits);

const PFLayer::Layer layer = rh.layer();
int iring = 0;
if (layer==PFLayer::HCAL_BARREL2 && abs(rh.positionREP().Eta())>0.34) iring= 1;

double rh_thresh = parameter( THRESH, layer, 0, iring );
if( rh_thresh <= 0.0 ) rh_thresh = 1.0; // in case of zero threshold do NOT normalize rechit energies

// int layer = rh.layer();

dist.clear();
Expand Down Expand Up @@ -1189,7 +1196,7 @@ PFClusterAlgo::buildPFClusters( const std::vector< unsigned >& topocluster,

// the current cell is the seed from the current photon.
if( rhindex == seedsintopocluster[ic] && seedexclusion ) {
frc = 1.;
frc = 1./rh_thresh;
#ifdef PFLOW_DEBUG
if(debug_) cout<<"this cell is a seed for the current photon"<<endl;
#endif
Expand All @@ -1203,7 +1210,12 @@ PFClusterAlgo::buildPFClusters( const std::vector< unsigned >& topocluster,
else {
// Compute the fractions of the cell energy to be assigned to
// each curpfclusters in the cluster.
frc = ener[ic] * exp ( - dist[ic]*dist[ic] / 2. );
// normalize the cluster energy to the rechit threshold
// to make clustering depend on *relative* energies of
// rechits to clusters
// this makes the numerical stability cut later not
// cut overly hard on the PS for instance
frc = (ener[ic]/rh_thresh) * exp ( - dist[ic]*dist[ic] / 2. );

#ifdef PFLOW_DEBUG
if(debug_) {
Expand All @@ -1228,7 +1240,7 @@ PFClusterAlgo::buildPFClusters( const std::vector< unsigned >& topocluster,
cout<<" frac["<<ic<<"] "<<frac[ic]<<" "<<fractot<<" "<<rh<<endl;
#endif

if( fractot )
if( fractot > 1e-20 || ( rh.detId() == curpfclusters[ic].seed() && fractot > 0.0 ) )
frac[ic] /= fractot;
else {
#ifdef PFLOW_DEBUG
Expand Down
8 changes: 8 additions & 0 deletions RecoParticleFlow/PFClusterProducer/test/BuildFile.xml
Expand Up @@ -6,3 +6,11 @@
<use name="FWCore/Utilities"/>
<flags EDM_PLUGIN="1"/>
</library>
<library name="PFClusterAnalyzer" file="PFClusterComparator.cc">
<use name="DataFormats/ParticleFlowReco"/>
<use name="FWCore/Framework"/>
<use name="FWCore/MessageLogger"/>
<use name="FWCore/ParameterSet"/>
<use name="FWCore/Utilities"/>
<flags EDM_PLUGIN="1"/>
</library>
184 changes: 184 additions & 0 deletions RecoParticleFlow/PFClusterProducer/test/PFClusterComparator.cc
@@ -0,0 +1,184 @@
#include "RecoParticleFlow/PFClusterProducer/test/PFClusterComparator.h"
#include "DataFormats/ParticleFlowReco/interface/PFCluster.h"

#include "FWCore/Framework/interface/ESHandle.h"

#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "FWCore/Utilities/interface/Exception.h"
#include "FWCore/Framework/interface/EventSetup.h"

#include <iomanip>

using namespace std;
using namespace edm;
using namespace reco;

PFClusterComparator::PFClusterComparator(const edm::ParameterSet& iConfig) {



inputTagPFClusters_
= iConfig.getParameter<InputTag>("PFClusters");

inputTagPFClustersCompare_
= iConfig.getParameter<InputTag>("PFClustersCompare");

verbose_ =
iConfig.getUntrackedParameter<bool>("verbose",false);

printBlocks_ =
iConfig.getUntrackedParameter<bool>("printBlocks",false);



LogDebug("PFClusterComparator")
<<" input collection : "<<inputTagPFClusters_ ;

}



PFClusterComparator::~PFClusterComparator() { }



void
PFClusterComparator::beginRun(const edm::Run& run,
const edm::EventSetup & es) { }


void PFClusterComparator::analyze(const Event& iEvent,
const EventSetup& iSetup) {

// get PFClusters

Handle<PFClusterCollection> pfClusters;
fetchCandidateCollection(pfClusters,
inputTagPFClusters_,
iEvent );

Handle<PFClusterCollection> pfClustersCompare;
fetchCandidateCollection(pfClustersCompare,
inputTagPFClustersCompare_,
iEvent );

// get PFClusters for isolation

std::cout << "There are " << pfClusters->size() << " PFClusters in the original cluster collection." << std::endl;
std::cout << "There are " << pfClustersCompare->size() << " PFClusters in the new cluster collection." << std::endl;

std::cout << std::flush << "---- COMPARING OLD TO NEW ----"
<< std::endl << std::flush;

for( unsigned i=0; i<pfClusters->size(); i++ ) {
const reco::PFCluster& cluster = pfClusters->at(i);
bool foundmatch = false;
for( unsigned k=0; k<pfClustersCompare->size(); ++k ) {
const reco::PFCluster& clustercomp = pfClustersCompare->at(k);
if( cluster.seed() == clustercomp.seed() ) {
foundmatch = true;
const double denergy = std::abs(cluster.energy() -
clustercomp.energy());
const double dcenergy = std::abs(cluster.correctedEnergy() -
clustercomp.correctedEnergy());
const double dx = std::abs(cluster.position().x() -
clustercomp.position().x());
const double dy = std::abs(cluster.position().y() -
clustercomp.position().y());
const double dz = std::abs(cluster.position().z() -
clustercomp.position().z());
if( denergy/std::abs(cluster.energy()) > 1e-5 ) {
std::cout << " " << cluster.seed()
<< " Energies different by larger than tolerance! "
<< "( "<< denergy << " )"
<< " Old: " << std::setprecision(7)
<< cluster.energy() << " GeV , New: "
<< clustercomp.energy() << " GeV" << std::endl;
}
if( dcenergy/std::abs(cluster.correctedEnergy()) > 1e-5 ) {
std::cout << " " << cluster.seed()
<< " Corrected energies different by larger than tolerance! "
<< "( "<< denergy << " )"
<< " Old: " << std::setprecision(7)
<< cluster.correctedEnergy() << " GeV , New: "
<< clustercomp.correctedEnergy() << " GeV" << std::endl;
}
std::cout << std::flush;
if( dx/std::abs(cluster.position().x()) > 1e-5 ) {
std::cout << "***" << cluster.seed()
<< " X's different by larger than tolerance! "
<< "( "<< dx << " )"
<< " Old: " << std::setprecision(7)
<< cluster.position().x() << " , New: "
<< clustercomp.position().x() << std::endl;
}
std::cout << std::flush;
if( dy/std::abs(cluster.position().y()) > 1e-5 ) {
std::cout << "---" << cluster.seed()
<< " Y's different by larger than tolerance! "
<< "( "<< dy << " )"
<< " Old: " << std::setprecision(7)
<< cluster.position().y() << " , New: "
<< clustercomp.position().y() << std::endl;
}
std::cout << std::flush;
if( dz/std::abs(cluster.position().z()) > 1e-5 ) {
std::cout << "+++" << cluster.seed()
<< " Z's different by larger than tolerance! "
<< "( "<< dz << " )"
<< " Old: " << std::setprecision(7)
<< cluster.position().z() << " , New: "
<< clustercomp.position().z() << std::endl;
}
std::cout << std::flush;
}
}
if( !foundmatch ) {
std::cout << "Seed in old clusters and not new: "
<< cluster << std::endl;
}
}

std::cout << std::flush << "---- COMPARING NEW TO OLD ----"
<< std::endl << std::flush;

for( unsigned i=0; i<pfClustersCompare->size(); i++ ) {
const reco::PFCluster& cluster = pfClustersCompare->at(i);
bool foundmatch = false;
for( unsigned k=0; k<pfClusters->size(); ++k ) {
const reco::PFCluster& clustercomp = pfClusters->at(k);
if( cluster.seed() == clustercomp.seed() ) {
foundmatch = true;

}
}
if( !foundmatch ) {
std::cout << "Seed in new clusters and not old: "
<< cluster << std::endl;
}
}
std::cout << std::flush;
}



void
PFClusterComparator::fetchCandidateCollection(Handle<reco::PFClusterCollection>& c,
const InputTag& tag,
const Event& iEvent) const {

bool found = iEvent.getByLabel(tag, c);

if(!found ) {
ostringstream err;
err<<" cannot get PFClusters: "
<<tag<<endl;
LogError("PFClusters")<<err.str();
throw cms::Exception( "MissingProduct", err.str());
}

}



DEFINE_FWK_MODULE(PFClusterComparator);
61 changes: 61 additions & 0 deletions RecoParticleFlow/PFClusterProducer/test/PFClusterComparator.h
@@ -0,0 +1,61 @@
#ifndef RecoParticleFlow_PFClusterProducer_PFClusterComparator_
#define RecoParticleFlow_PFClusterProducer_PFClusterComparator_

// system include files
#include <memory>
#include <string>
#include <iostream>

// user include files
#include "FWCore/Framework/interface/Frameworkfwd.h"
#include "FWCore/Framework/interface/EDAnalyzer.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"

#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/MakerMacros.h"

#include "DataFormats/ParticleFlowReco/interface/PFClusterFwd.h"

/**\class PFClusterAnalyzer
\brief test analyzer for PFClusters
*/




class PFClusterComparator : public edm::EDAnalyzer {
public:

explicit PFClusterComparator(const edm::ParameterSet&);

~PFClusterComparator();

virtual void analyze(const edm::Event&, const edm::EventSetup&);

virtual void beginRun(const edm::Run & r, const edm::EventSetup & c);

private:

void
fetchCandidateCollection(edm::Handle<reco::PFClusterCollection>& c,
const edm::InputTag& tag,
const edm::Event& iSetup) const;

/* void printElementsInBlocks(const reco::PFCluster& cluster, */
/* std::ostream& out=std::cout) const; */



/// PFClusters in which we'll look for pile up particles
edm::InputTag inputTagPFClusters_;
edm::InputTag inputTagPFClustersCompare_;

/// verbose ?
bool verbose_;

/// print the blocks associated to a given candidate ?
bool printBlocks_;

};

#endif

0 comments on commit 4d93095

Please sign in to comment.