13 #include "Math/GenVector/VectorUtil.h" 35 typedef std::pair<reco::CaloClusterPtr::key_type,reco::CaloClusterPtr>
EEPSPair;
36 typedef std::binary_function<
const CalibClusterPtr&,
37 const CalibClusterPtr&,
38 bool> ClusBinaryFunction;
40 typedef std::unary_function<
const CalibClusterPtr&,
41 bool> ClusUnaryFunction;
43 bool sortByKey(
const EEPSPair&
a,
const EEPSPair&
b) {
44 return a.first < b.first;
47 inline double getPFClusterEnergy(
const PFClusterPtr&
p) {
51 inline double ptFast(
const double energy,
54 const auto v = position - origin;
58 struct SumPSEnergy :
public std::binary_function<double,
63 double operator()(
double a,
64 const PFClusterPtr&
b) {
65 return a + (_thelayer == b->layer())*b->energy();
69 struct GreaterByE :
public ClusBinaryFunction {
70 bool operator()(
const CalibClusterPtr& x,
71 const CalibClusterPtr& y) {
72 return x->energy() > y->energy() ;
76 struct GreaterByEt :
public ClusBinaryFunction {
77 bool operator()(
const CalibClusterPtr& x,
78 const CalibClusterPtr& y) {
80 const double xpt =
ptFast(x->energy(),x->the_ptr()->position(),zero);
81 const double ypt =
ptFast(y->energy(),y->the_ptr()->position(),zero);
86 struct IsASeed :
public ClusUnaryFunction {
89 IsASeed(
double thresh,
bool useETcut =
false) :
90 threshold(thresh), cutET(useETcut) {}
91 bool operator()(
const CalibClusterPtr& x) {
93 double e_or_et = x->energy();
94 if( cutET ) e_or_et =
ptFast(e_or_et,x->the_ptr()->position(),zero);
99 struct IsLinkedByRecHit :
public ClusUnaryFunction {
100 const CalibClusterPtr the_seed;
101 const double _threshold, _majority;
102 const double _maxSatelliteDEta, _maxSatelliteDPhi;
103 double x_rechits_tot=0;
104 double x_rechits_match=0;
105 IsLinkedByRecHit(
const CalibClusterPtr&
s,
const double threshold,
106 const double majority,
const double maxDEta,
108 the_seed(s),_threshold(threshold),_majority(majority),
109 _maxSatelliteDEta(maxDEta), _maxSatelliteDPhi(maxDPhi) {}
110 bool operator()(
const CalibClusterPtr& x) {
111 if( the_seed->energy_nocalib() < _threshold )
return false;
112 const double dEta =
std::abs(the_seed->eta()-x->eta());
114 std::abs(TVector2::Phi_mpi_pi(the_seed->phi() - x->phi()));
115 if( _maxSatelliteDEta < dEta || _maxSatelliteDPhi < dPhi)
return false;
117 const auto& seedHitsAndFractions =
118 the_seed->the_ptr()->hitsAndFractions();
119 const auto& xHitsAndFractions =
120 x->the_ptr()->hitsAndFractions();
121 x_rechits_tot = xHitsAndFractions.size();
122 x_rechits_match = 0.0;
123 for(
const std::pair<DetId, float>& seedHit : seedHitsAndFractions ) {
124 for(
const std::pair<DetId, float>& xHit : xHitsAndFractions ) {
125 if( seedHit.first == xHit.first ) {
126 x_rechits_match += 1.0;
130 return x_rechits_match/x_rechits_tot > _majority;
134 struct IsClustered :
public ClusUnaryFunction {
135 const CalibClusterPtr the_seed;
138 double etawidthSuperCluster_ = .0 , phiwidthSuperCluster_ = .0;
139 IsClustered(
const CalibClusterPtr
s,
141 const bool dyn_dphi) :
142 the_seed(s), _type(ct), dynamic_dphi(dyn_dphi) {}
143 bool operator()(
const CalibClusterPtr& x) {
145 std::abs(TVector2::Phi_mpi_pi(the_seed->phi() - x->phi()));
146 const bool passes_dphi =
147 ( (!dynamic_dphi && dphi < phiwidthSuperCluster_ ) ||
156 return (
std::abs(the_seed->eta()-x->eta())<etawidthSuperCluster_ &&
160 return ( passes_dphi &&
195 regr_->setTokens(regconf, cc);
209 regr_->setEventSetup(setup);
245 regr_->setEvent(iEvent);
256 for (
size_t i = 0;
i < clusters.size(); ++
i ){
257 auto cluster = clusters.ptrAt(
i);
259 <<
"Loading PFCluster i="<<cluster.key()
260 <<
" energy="<<cluster->energy()<<std::endl;
263 if( cluster->caloID().detectors() == 0 &&
264 cluster->hitsAndFractions().empty() )
continue;
267 switch( cluster->layer() ) {
294 if (!barrelRecHitsHandle.
isValid()) {
296 <<
"If you use OOT photons, need to specify proper barrel rec hit collection";
302 if (!endcapRecHitsHandle.
isValid()) {
304 <<
"If you use OOT photons, need to specify proper endcap rec hit collection";
322 std::stable_partition(clusters.begin(),clusters.end(),seedable);
327 while( std::any_of(clusters.cbegin(), clusters.cend(), seedable) ) {
340 switch( seed->the_ptr()->layer() ) {
346 <<
" in the ECAL barrel!";
355 <<
" in the ECAL endcap!" << std::endl;
367 auto not_clustered = std::stable_partition(clusters.begin(),clusters.end(),
368 IsClusteredWithSeed);
372 not_clustered = std::stable_partition(not_clustered,clusters.end(),
373 MatchesSeedByRecHit);
377 edm::LogInfo(
"PFClustering") <<
"Dumping cluster detail";
379 <<
"\tPassed seed: e = " << seed->energy_nocalib()
380 <<
" eta = " << seed->eta() <<
" phi = " << seed->phi()
382 for(
auto clus = clusters.cbegin(); clus != not_clustered; ++clus ) {
384 <<
"\t\tClustered cluster: e = " << (*clus)->energy_nocalib()
385 <<
" eta = " << (*clus)->eta() <<
" phi = " << (*clus)->phi()
388 for(
auto clus = not_clustered; clus != clusters.end(); ++clus ) {
390 <<
"\tNon-Clustered cluster: e = " << (*clus)->energy_nocalib()
391 <<
" eta = " << (*clus)->eta() <<
" phi = " << (*clus)->phi()
396 if (not_clustered == clusters.begin()) {
398 clusters.erase(clusters.begin());
403 <<
"Cluster is not seedable!" << std::endl
404 <<
"\tNon-Clustered cluster: e = " << (*not_clustered)->energy_nocalib()
405 <<
" eta = " << (*not_clustered)->eta() <<
" phi = " << (*not_clustered)->phi()
413 clusters.erase(clusters.begin(),not_clustered);
415 std::vector<const reco::PFCluster*> bare_ptrs;
418 corrSCEnergy(0), corrPS1Energy(0), corrPS2Energy(0),
419 ePS1(0), ePS2(0), energyweight(0), energyweighttot(0);
420 std::vector<double> ps1_energies, ps2_energies;
421 int condP1(1), condP2(1);
422 for(
auto& clus : clustered ) {
424 energyweight = clus->energy_nocalib();
425 bare_ptrs.push_back(clus->the_ptr().get());
430 ps1_energies.clear();
431 ps2_energies.clear();
434 const auto clustops = std::equal_range(
EEtoPS_->begin(),
438 for(
auto i_ps = clustops.first; i_ps != clustops.second; ++i_ps) {
443 switch( psclus->
layer() ) {
445 ps1_energies.push_back(psclus->
energy());
446 for (
auto const& recH : recH_Frac){
447 ESDetId strip1 = recH.recHitRef()->detId();
452 if(status_p1->getStatusCode() == 0) condP1 = 0;
457 ps2_energies.push_back(psclus->
energy());
458 for (
auto const& recH : recH_Frac){
459 ESDetId strip2 = recH.recHitRef()->detId();
462 if(status_p2->getStatusCode() == 0) condP2 = 0;
470 if(condP1 == 1) ePS1 = -1.;
471 if(condP2 == 1) ePS2 = -1.;
473 ps1_energies,ps2_energies,
478 if(ePS1 == -1.) ePS1 = 0;
479 if(ePS2 == -1.) ePS2 = 0;
485 energyweight = clus->energy() - ePS1 - ePS2;
488 energyweight = clus->energy();
494 posX += energyweight * cluspos.X();
495 posY += energyweight * cluspos.Y();
496 posZ += energyweight * cluspos.Z();
498 energyweighttot += energyweight;
499 corrSCEnergy += clus->energy();
500 corrPS1Energy += ePS1;
501 corrPS2Energy += ePS2;
503 posX /= energyweighttot;
504 posY /= energyweighttot;
505 posZ /= energyweighttot;
510 new_sc.
setSeed(clustered.front()->the_ptr());
514 for(
const auto& clus : clustered ) {
517 auto& hits_and_fractions = clus->the_ptr()->hitsAndFractions();
518 for(
auto& hit_and_fraction : hits_and_fractions ) {
524 const auto clustops = std::equal_range(
EEtoPS_->begin(),
531 for(
auto i_ps = clustops.first; i_ps != clustops.second; ++i_ps) {
543 <<
"Found a PS cluster matched to more than one EE cluster!" 544 << std::endl << std::hex << psclus.
get() <<
" == " 545 << found_pscluster->get() <<
std::dec << std::endl;
562 regr_->modifyObject(new_sc);
574 switch( seed->the_ptr()->layer() ) {
PFLayer::Layer layer() const
cluster layer, see PFLayer.h in this directory
T getParameter(std::string const &) const
const math::XYZPoint & position() const
cluster centroid position
CaloCluster_iterator preshowerClustersBegin() const
fist iterator over PreshowerCluster constituents
double etawidthSuperClusterBarrel_
PFECALSuperClusterAlgo()
constructor
const ESChannelStatus * channelStatus_
double phiwidthSuperClusterEndcap_
void addHitAndFraction(DetId id, float fraction)
edm::EDGetTokenT< edm::View< reco::PFCluster > > inputTagPFClusters_
std::unique_ptr< SCEnergyCorrectorSemiParm > regr_
bool inDynamicDPhiWindow(const float seedEta, const float seedPhi, const float ClustE, const float ClusEta, const float clusPhi)
bool getByToken(EDGetToken token, Handle< PROD > &result) const
std::unique_ptr< reco::SuperClusterCollection > superClustersEB_
void setPreshowerEnergyPlane2(double preshowerEnergy2)
edm::EDGetTokenT< EcalRecHitCollection > inputTagBarrelRecHits_
const self & getMap() const
T const * get() const
Returns C++ pointer to the item.
double threshPFClusterSeedBarrel_
bool applyCrackCorrections_
double pflowPhiWidth() const
std::vector< EcalRecHit >::const_iterator const_iterator
def setup(process, global_tag, zero_tesla=False)
const reco::BeamSpot * beamSpot_
const reco::PFCluster::EEtoPSAssociation * EEtoPS_
void setSeed(const CaloClusterPtr &r)
list of used xtals by DetId // now inherited by CaloCluster
void find(edm::Handle< EcalRecHitCollection > &hits, DetId thisDet, std::vector< EcalRecHitCollection::const_iterator > &hit, bool debug=false)
void setPhiWidth(double pw)
double etawidthSuperClusterEndcap_
double threshPFClusterSeedEndcap_
double pflowEtaWidth() const
double ptFast(const double energy, const math::XYZPoint &position, const math::XYZPoint &origin)
std::unique_ptr< reco::SuperClusterCollection > superClustersEE_
edm::EDGetTokenT< EcalRecHitCollection > inputTagEndcapRecHits_
void loadAndSortPFClusters(const edm::Event &evt)
void update(const edm::EventSetup &)
void setEtaWidth(double ew)
std::vector< CalibratedClusterPtr > CalibratedClusterPtrVector
std::shared_ptr< CalibratedPFCluster > CalibratedClusterPtr
MVATrainerComputer * calib
void buildAllSuperClusters(CalibratedClusterPtrVector &, double seedthresh)
std::shared_ptr< PFEnergyCalibration > _pfEnergyCalibration
void setTokens(const edm::ParameterSet &, edm::ConsumesCollector &&)
bool doSatelliteClusterMerge_
void setCorrectedEnergy(double cenergy)
const_iterator find(uint32_t rawId) const
std::vector< std::pair< CaloClusterPtr::key_type, edm::Ptr< PFCluster > > > EEtoPSAssociation
Abs< T >::type abs(const T &t)
CalibratedClusterPtrVector _clustersEE
double threshSuperClusterEt_
double energy() const
cluster energy
double satelliteThreshold_
double fractionForMajority_
double threshPFClusterBarrel_
double energy() const
cluster energy
double phiwidthSuperClusterBarrel_
double rawEnergy() const
raw uncorrected energy (sum of energies of component BasicClusters)
edm::EDGetTokenT< reco::PFCluster::EEtoPSAssociation > inputTagPFClustersES_
T const * product() const
edm::EDGetTokenT< reco::BeamSpot > inputTagBeamSpot_
XYZPointD XYZPoint
point in space with cartesian internal representation
reco::PFCluster::EEtoPSAssociation::value_type EEPSPair
void addPreshowerCluster(const CaloClusterPtr &r)
add reference to constituent BasicCluster
bool sortByKey(const EEPSPair &a, const EEPSPair &b)
clustering_type _clustype
double threshPFClusterEndcap_
void buildSuperCluster(CalibratedClusterPtr &, CalibratedClusterPtrVector &)
std::vector< Item >::const_iterator const_iterator
bool GreaterByE(const T &a1, const T &a2)
void addCluster(const CaloClusterPtr &r)
add reference to constituent BasicCluster
iterator find(key_type k)
static int position[264][3]
const std::vector< reco::PFRecHitFraction > & recHitFractions() const
vector of rechit fractions
const Point & position() const
position
const CaloClusterPtr & seed() const
seed BasicCluster
const EcalRecHitCollection * barrelRecHits_
T const * product() const
void setPreshowerEnergyPlane1(double preshowerEnergy1)
void setPFClusterCalibration(const std::shared_ptr< PFEnergyCalibration > &)
CalibratedClusterPtrVector _clustersEB
const EcalRecHitCollection * endcapRecHits_
bool inMustache(const float maxEta, const float maxPhi, const float ClustE, const float ClusEta, const float ClusPhi)
CaloCluster_iterator preshowerClustersEnd() const
last iterator over PreshowerCluster constituents
void setPreshowerEnergy(double preshowerEnergy)