104 auto const & hits = *
input;
105 std::vector<unsigned > ordered_hits(hits.size());
106 for(
unsigned i = 0;
i < hits.size(); ++
i ) ordered_hits[
i]=
i;
107 std::sort(ordered_hits.begin(),ordered_hits.end(),[&](
unsigned i,
unsigned j) {
return hits[
i].energy()>hits[
j].energy();});
109 for(
const auto&
idx : ordered_hits ) {
110 const unsigned i =
idx;
111 if( !mask[i] )
continue;
113 int hitlayer = (int)rechit.
layer();
119 if( rechit.
energy() < clean._singleSpikeThresh )
continue;
120 const double rhenergy = rechit.
energy();
123 double surroundingEnergy = rechit.
energy();
124 double neighbourEnergy = 0.0;
125 double layerEnergy = 0.0;
126 for(
auto k : neighbours4 ) {
127 if( !mask[
k] )
continue;
128 const auto & neighbour = hits[
k];
129 const double sum = neighbour.energy();
130 surroundingEnergy += sum;
131 neighbourEnergy += sum;
132 layerEnergy += neighbour.energy();
137 const double fraction1 = surroundingEnergy/rhenergy;
140 const double f1Cut = ( clean._minS4S1_a*std::log10(rechit.
energy()) +
142 if( fraction1 < f1Cut ) {
146 std::pair<double,double> dcr = dCrack(phi,eta);
151 ( (aeta < 2.85 && dcrmin > 1.0) ||
152 (rhenergy > clean._eneThreshMod*clean._singleSpikeThresh &&
153 fraction1 < f1Cut/clean._fracThreshMod ) ) ) {
158 if( mask[i] && rhenergy > clean._doubleSpikeThresh ) {
160 double surroundingEnergyi = 0.0;
161 double enmax = -999.0;
162 unsigned int mostEnergeticNeighbour=0;
164 for(
auto k : neighbours4i ) {
165 if( !mask[k] )
continue;
166 const auto & neighbour = hits[
k];
167 const double nenergy = neighbour.energy();
168 surroundingEnergyi += nenergy;
169 if( nenergy > enmax ) {
171 mostEnergeticNeighbour =
k;
176 double surroundingEnergyj = 0.0;
177 auto const & neighbours4j =
178 hits[mostEnergeticNeighbour].neighbours4();
179 for(
auto k : neighbours4j ) {
181 surroundingEnergyj += hits[
k].energy();
184 const double surroundingEnergyFraction =
185 (surroundingEnergyi+surroundingEnergyj) /
186 (rechit.
energy()+hits[mostEnergeticNeighbour].energy()) - 1.;
187 if ( surroundingEnergyFraction < clean._doubleSpikeS6S2 ) {
191 std::pair<double,double> dcr = dCrack(phi,eta);
196 ( (aeta < 2.85 && dcrmin > 1.0) ||
197 (rhenergy > clean._eneThreshMod*clean._doubleSpikeThresh &&
198 surroundingEnergyFraction < clean._doubleSpikeS6S2/clean._fracThreshMod )
201 mask[mostEnergeticNeighbour] =
false;
float phi() const
momentum azimuthal angle
void clean(const edm::Handle< reco::PFRecHitCollection > &input, std::vector< bool > &mask)
static std::string const input
PFLayer::Layer layer() const
rechit layer
U second(std::pair< T, U > const &p)
Particle flow rechit (rechit + geometry and topology information). See clustering algorithm in PFClus...
Abs< T >::type abs(const T &t)
float energy() const
rechit energy
tuple idx
DEBUGGING if hasattr(process,"trackMonIterativeTracking2012"): print "trackMonIterativeTracking2012 D...
std::unordered_map< int, spike_cleaning > _thresholds
float eta() const
momentum pseudorapidity
Neighbours neighbours4() const
RhoEtaPhi const & positionREP() const