Skip to content

Commit 4d93095

Browse files
committed
Merge pull request cms-sw#2490 from lgray/pfclusteralgo_stability
Reco updates -- Protection to stabilize the behavior of PFClusterAlgo
2 parents 1910581 + f88a2d9 commit 4d93095

File tree

5 files changed

+269
-4
lines changed

5 files changed

+269
-4
lines changed

RecoParticleFlow/PFClusterProducer/python/particleFlowClusterPS_cfi.py

+1-1
Original file line numberDiff line numberDiff line change
@@ -35,7 +35,7 @@
3535
# n neighbours in PS
3636
nNeighbours = cms.int32(8),
3737
# sigma of the shower in preshower
38-
showerSigma = cms.double(0.2),
38+
showerSigma = cms.double(0.3),
3939
# n crystals for position calculation in PS
4040
posCalcNCrystal = cms.int32(-1),
4141
# use cells with common corner to build topo-clusters

RecoParticleFlow/PFClusterProducer/src/PFClusterAlgo.cc

+15-3
Original file line numberDiff line numberDiff line change
@@ -1122,6 +1122,13 @@ PFClusterAlgo::buildPFClusters( const std::vector< unsigned >& topocluster,
11221122

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

1125+
const PFLayer::Layer layer = rh.layer();
1126+
int iring = 0;
1127+
if (layer==PFLayer::HCAL_BARREL2 && abs(rh.positionREP().Eta())>0.34) iring= 1;
1128+
1129+
double rh_thresh = parameter( THRESH, layer, 0, iring );
1130+
if( rh_thresh <= 0.0 ) rh_thresh = 1.0; // in case of zero threshold do NOT normalize rechit energies
1131+
11251132
// int layer = rh.layer();
11261133

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

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

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

1231-
if( fractot )
1243+
if( fractot > 1e-20 || ( rh.detId() == curpfclusters[ic].seed() && fractot > 0.0 ) )
12321244
frac[ic] /= fractot;
12331245
else {
12341246
#ifdef PFLOW_DEBUG

RecoParticleFlow/PFClusterProducer/test/BuildFile.xml

+8
Original file line numberDiff line numberDiff line change
@@ -6,3 +6,11 @@
66
<use name="FWCore/Utilities"/>
77
<flags EDM_PLUGIN="1"/>
88
</library>
9+
<library name="PFClusterAnalyzer" file="PFClusterComparator.cc">
10+
<use name="DataFormats/ParticleFlowReco"/>
11+
<use name="FWCore/Framework"/>
12+
<use name="FWCore/MessageLogger"/>
13+
<use name="FWCore/ParameterSet"/>
14+
<use name="FWCore/Utilities"/>
15+
<flags EDM_PLUGIN="1"/>
16+
</library>
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,184 @@
1+
#include "RecoParticleFlow/PFClusterProducer/test/PFClusterComparator.h"
2+
#include "DataFormats/ParticleFlowReco/interface/PFCluster.h"
3+
4+
#include "FWCore/Framework/interface/ESHandle.h"
5+
6+
#include "FWCore/MessageLogger/interface/MessageLogger.h"
7+
#include "FWCore/Utilities/interface/Exception.h"
8+
#include "FWCore/Framework/interface/EventSetup.h"
9+
10+
#include <iomanip>
11+
12+
using namespace std;
13+
using namespace edm;
14+
using namespace reco;
15+
16+
PFClusterComparator::PFClusterComparator(const edm::ParameterSet& iConfig) {
17+
18+
19+
20+
inputTagPFClusters_
21+
= iConfig.getParameter<InputTag>("PFClusters");
22+
23+
inputTagPFClustersCompare_
24+
= iConfig.getParameter<InputTag>("PFClustersCompare");
25+
26+
verbose_ =
27+
iConfig.getUntrackedParameter<bool>("verbose",false);
28+
29+
printBlocks_ =
30+
iConfig.getUntrackedParameter<bool>("printBlocks",false);
31+
32+
33+
34+
LogDebug("PFClusterComparator")
35+
<<" input collection : "<<inputTagPFClusters_ ;
36+
37+
}
38+
39+
40+
41+
PFClusterComparator::~PFClusterComparator() { }
42+
43+
44+
45+
void
46+
PFClusterComparator::beginRun(const edm::Run& run,
47+
const edm::EventSetup & es) { }
48+
49+
50+
void PFClusterComparator::analyze(const Event& iEvent,
51+
const EventSetup& iSetup) {
52+
53+
// get PFClusters
54+
55+
Handle<PFClusterCollection> pfClusters;
56+
fetchCandidateCollection(pfClusters,
57+
inputTagPFClusters_,
58+
iEvent );
59+
60+
Handle<PFClusterCollection> pfClustersCompare;
61+
fetchCandidateCollection(pfClustersCompare,
62+
inputTagPFClustersCompare_,
63+
iEvent );
64+
65+
// get PFClusters for isolation
66+
67+
std::cout << "There are " << pfClusters->size() << " PFClusters in the original cluster collection." << std::endl;
68+
std::cout << "There are " << pfClustersCompare->size() << " PFClusters in the new cluster collection." << std::endl;
69+
70+
std::cout << std::flush << "---- COMPARING OLD TO NEW ----"
71+
<< std::endl << std::flush;
72+
73+
for( unsigned i=0; i<pfClusters->size(); i++ ) {
74+
const reco::PFCluster& cluster = pfClusters->at(i);
75+
bool foundmatch = false;
76+
for( unsigned k=0; k<pfClustersCompare->size(); ++k ) {
77+
const reco::PFCluster& clustercomp = pfClustersCompare->at(k);
78+
if( cluster.seed() == clustercomp.seed() ) {
79+
foundmatch = true;
80+
const double denergy = std::abs(cluster.energy() -
81+
clustercomp.energy());
82+
const double dcenergy = std::abs(cluster.correctedEnergy() -
83+
clustercomp.correctedEnergy());
84+
const double dx = std::abs(cluster.position().x() -
85+
clustercomp.position().x());
86+
const double dy = std::abs(cluster.position().y() -
87+
clustercomp.position().y());
88+
const double dz = std::abs(cluster.position().z() -
89+
clustercomp.position().z());
90+
if( denergy/std::abs(cluster.energy()) > 1e-5 ) {
91+
std::cout << " " << cluster.seed()
92+
<< " Energies different by larger than tolerance! "
93+
<< "( "<< denergy << " )"
94+
<< " Old: " << std::setprecision(7)
95+
<< cluster.energy() << " GeV , New: "
96+
<< clustercomp.energy() << " GeV" << std::endl;
97+
}
98+
if( dcenergy/std::abs(cluster.correctedEnergy()) > 1e-5 ) {
99+
std::cout << " " << cluster.seed()
100+
<< " Corrected energies different by larger than tolerance! "
101+
<< "( "<< denergy << " )"
102+
<< " Old: " << std::setprecision(7)
103+
<< cluster.correctedEnergy() << " GeV , New: "
104+
<< clustercomp.correctedEnergy() << " GeV" << std::endl;
105+
}
106+
std::cout << std::flush;
107+
if( dx/std::abs(cluster.position().x()) > 1e-5 ) {
108+
std::cout << "***" << cluster.seed()
109+
<< " X's different by larger than tolerance! "
110+
<< "( "<< dx << " )"
111+
<< " Old: " << std::setprecision(7)
112+
<< cluster.position().x() << " , New: "
113+
<< clustercomp.position().x() << std::endl;
114+
}
115+
std::cout << std::flush;
116+
if( dy/std::abs(cluster.position().y()) > 1e-5 ) {
117+
std::cout << "---" << cluster.seed()
118+
<< " Y's different by larger than tolerance! "
119+
<< "( "<< dy << " )"
120+
<< " Old: " << std::setprecision(7)
121+
<< cluster.position().y() << " , New: "
122+
<< clustercomp.position().y() << std::endl;
123+
}
124+
std::cout << std::flush;
125+
if( dz/std::abs(cluster.position().z()) > 1e-5 ) {
126+
std::cout << "+++" << cluster.seed()
127+
<< " Z's different by larger than tolerance! "
128+
<< "( "<< dz << " )"
129+
<< " Old: " << std::setprecision(7)
130+
<< cluster.position().z() << " , New: "
131+
<< clustercomp.position().z() << std::endl;
132+
}
133+
std::cout << std::flush;
134+
}
135+
}
136+
if( !foundmatch ) {
137+
std::cout << "Seed in old clusters and not new: "
138+
<< cluster << std::endl;
139+
}
140+
}
141+
142+
std::cout << std::flush << "---- COMPARING NEW TO OLD ----"
143+
<< std::endl << std::flush;
144+
145+
for( unsigned i=0; i<pfClustersCompare->size(); i++ ) {
146+
const reco::PFCluster& cluster = pfClustersCompare->at(i);
147+
bool foundmatch = false;
148+
for( unsigned k=0; k<pfClusters->size(); ++k ) {
149+
const reco::PFCluster& clustercomp = pfClusters->at(k);
150+
if( cluster.seed() == clustercomp.seed() ) {
151+
foundmatch = true;
152+
153+
}
154+
}
155+
if( !foundmatch ) {
156+
std::cout << "Seed in new clusters and not old: "
157+
<< cluster << std::endl;
158+
}
159+
}
160+
std::cout << std::flush;
161+
}
162+
163+
164+
165+
void
166+
PFClusterComparator::fetchCandidateCollection(Handle<reco::PFClusterCollection>& c,
167+
const InputTag& tag,
168+
const Event& iEvent) const {
169+
170+
bool found = iEvent.getByLabel(tag, c);
171+
172+
if(!found ) {
173+
ostringstream err;
174+
err<<" cannot get PFClusters: "
175+
<<tag<<endl;
176+
LogError("PFClusters")<<err.str();
177+
throw cms::Exception( "MissingProduct", err.str());
178+
}
179+
180+
}
181+
182+
183+
184+
DEFINE_FWK_MODULE(PFClusterComparator);
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,61 @@
1+
#ifndef RecoParticleFlow_PFClusterProducer_PFClusterComparator_
2+
#define RecoParticleFlow_PFClusterProducer_PFClusterComparator_
3+
4+
// system include files
5+
#include <memory>
6+
#include <string>
7+
#include <iostream>
8+
9+
// user include files
10+
#include "FWCore/Framework/interface/Frameworkfwd.h"
11+
#include "FWCore/Framework/interface/EDAnalyzer.h"
12+
#include "FWCore/ParameterSet/interface/ParameterSet.h"
13+
14+
#include "FWCore/Framework/interface/Event.h"
15+
#include "FWCore/Framework/interface/MakerMacros.h"
16+
17+
#include "DataFormats/ParticleFlowReco/interface/PFClusterFwd.h"
18+
19+
/**\class PFClusterAnalyzer
20+
\brief test analyzer for PFClusters
21+
*/
22+
23+
24+
25+
26+
class PFClusterComparator : public edm::EDAnalyzer {
27+
public:
28+
29+
explicit PFClusterComparator(const edm::ParameterSet&);
30+
31+
~PFClusterComparator();
32+
33+
virtual void analyze(const edm::Event&, const edm::EventSetup&);
34+
35+
virtual void beginRun(const edm::Run & r, const edm::EventSetup & c);
36+
37+
private:
38+
39+
void
40+
fetchCandidateCollection(edm::Handle<reco::PFClusterCollection>& c,
41+
const edm::InputTag& tag,
42+
const edm::Event& iSetup) const;
43+
44+
/* void printElementsInBlocks(const reco::PFCluster& cluster, */
45+
/* std::ostream& out=std::cout) const; */
46+
47+
48+
49+
/// PFClusters in which we'll look for pile up particles
50+
edm::InputTag inputTagPFClusters_;
51+
edm::InputTag inputTagPFClustersCompare_;
52+
53+
/// verbose ?
54+
bool verbose_;
55+
56+
/// print the blocks associated to a given candidate ?
57+
bool printBlocks_;
58+
59+
};
60+
61+
#endif

0 commit comments

Comments
 (0)