CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
List of all members | Public Types | Public Member Functions | Private Member Functions | Private Attributes | Static Private Attributes | Friends
PFClusterAlgo Class Reference

#include <PFClusterAlgo.h>

Public Types

enum  Color { NONE =0, SEED, SPECIAL }
 
typedef std::multimap< double,
unsigned >::iterator 
EH
 
typedef std::map< unsigned,
unsigned >::const_iterator 
IDH
 
enum  Parameter {
  THRESH, SEED_THRESH, PT_THRESH, SEED_PT_THRESH,
  CLEAN_THRESH, CLEAN_S4S1, DOUBLESPIKE_THRESH, DOUBLESPIKE_S6S2
}
 
typedef edm::Handle< edm::View
< reco::PFCluster > > 
PFClusterHandle
 
typedef edm::Handle
< reco::PFRecHitCollection
PFRecHitHandle
 
typedef
edm::StrictWeakOrdering
< reco::PFRecHit
PFStrictWeakOrdering
 
enum  PositionCalcType { EGPositionCalc, EGPositionFormula, PFPositionCalc, kNotDefined }
 
enum  SeedState { UNKNOWN =-1, NO =0, YES =1, CLEAN =2 }
 
typedef edm::SortedCollection
< reco::PFRecHit
SortedPFRecHitCollection
 

Public Member Functions

std::auto_ptr< std::vector
< reco::PFCluster > > & 
clusters ()
 
unsigned color (unsigned rhi) const
 
void doClustering (const reco::PFRecHitCollection &rechits)
 perform clustering More...
 
void doClustering (const reco::PFRecHitCollection &rechits, const std::vector< bool > &mask)
 
void doClustering (const PFRecHitHandle &rechitsHandle)
 perform clustering in full framework More...
 
void doClustering (const PFRecHitHandle &rechitsHandle, const std::vector< bool > &mask)
 
void enableDebugging (bool debug)
 set hits on which clustering will be performed More...
 
bool isSeed (unsigned rhi) const
 
bool masked (unsigned rhi) const
 
int nNeighbours () const
 get number of neighbours for More...
 
double parameter (Parameter paramtype, PFLayer::Layer layer, unsigned iCoeff=0, int iring0=0) const
 
 PFClusterAlgo ()
 constructor More...
 
int posCalcNCrystal () const
 get number of crystals for position calculation (-1 all,5, or 9) More...
 
double posCalcP1 () const
 get p1 for position calculation More...
 
const reco::PFRecHitrechit (unsigned i, const reco::PFRecHitCollection &rechits)
 
More...
 
std::auto_ptr< std::vector
< reco::PFRecHit > > & 
rechitsCleaned ()
 
void setCleanRBXandHPDs (bool cleanRBXandHPDs)
 Activate cleaning of HCAL RBX's and HPD's. More...
 
void setEBGeom (const CaloSubdetectorGeometry *esh)
 
void setEEGeom (const CaloSubdetectorGeometry *esh)
 
void setEGammaPosCalc (const edm::ParameterSet &conf)
 
void setHistos (TFile *file, TH2F *hB, TH2F *hE)
 set endcap clean threshold More...
 
void setNNeighbours (int n)
 set number of neighbours for More...
 
void setPosCalcNCrystal (int n)
 set number of crystals for position calculation (-1 all,5, or 9) More...
 
void setPosCalcP1 (double p1)
 set p1 for position calculation More...
 
void setPosCalcW0 (double w0)
 
void setPositionCalcType (const PositionCalcType &t)
 
void setPreshowerGeom (const CaloSubdetectorGeometry *esh)
 
void setS4S1CleanBarrel (const std::vector< double > &coeffs)
 
void setS4S1CleanEndcap (const std::vector< double > &coeffs)
 
void setS6S2DoubleSpikeBarrel (double cut)
 
void setS6S2DoubleSpikeEndcap (double cut)
 
void setShowerSigma (double sigma)
 set shower sigma for More...
 
void setThreshBarrel (double thresh)
 setters -------------------------------------------------—— More...
 
void setThreshCleanBarrel (double thresh)
 set barrel clean threshold More...
 
void setThreshCleanEndcap (double thresh)
 set endcap clean threshold More...
 
void setThreshDoubleSpikeBarrel (double thresh)
 set endcap thresholds for double spike cleaning More...
 
void setThreshDoubleSpikeEndcap (double thresh)
 set endcap thresholds for double spike cleaning More...
 
void setThreshEndcap (double thresh)
 set endcap threshold More...
 
void setThreshPtBarrel (double thresh)
 
void setThreshPtEndcap (double thresh)
 
void setThreshPtSeedBarrel (double thresh)
 
void setThreshPtSeedEndcap (double thresh)
 
void setThreshSeedBarrel (double thresh)
 set barrel seed threshold More...
 
void setThreshSeedEndcap (double thresh)
 set endcap seed threshold More...
 
void setUseCornerCells (bool usecornercells)
 activate use of cells with a common corner to build topo-clusters More...
 
double showerSigma () const
 get shower sigma More...
 
double threshBarrel () const
 getters -------------------------------------------------—— More...
 
double threshEndcap () const
 get endcap threshold More...
 
double threshSeedBarrel () const
 get barrel seed threshold More...
 
double threshSeedEndcap () const
 get endcap seed threshold More...
 
void write ()
 write histos More...
 
virtual ~PFClusterAlgo ()
 destructor More...
 

Private Member Functions

void buildPFClusters (const std::vector< unsigned > &cluster, const reco::PFRecHitCollection &rechits)
 build PFClusters from a topocluster More...
 
void buildTopoCluster (std::vector< unsigned > &cluster, unsigned rhi, const reco::PFRecHitCollection &rechits)
 build a topocluster (recursive) More...
 
void buildTopoClusters (const reco::PFRecHitCollection &rechits)
 build topoclusters around seeds More...
 
void calculateClusterPosition (reco::PFCluster &cluster, reco::PFCluster &clusterwodepthcor, bool depcor=true, int posCalcNCrystal=0)
 calculate position of a cluster More...
 
void cleanRBXAndHPD (const reco::PFRecHitCollection &rechits)
 Clean HCAL readout box noise and HPD discharge. More...
 
reco::PFRecHitRef createRecHitRef (const reco::PFRecHitCollection &rechits, unsigned rhi)
 
std::pair< double, double > dCrack (double phi, double eta)
 distance to a crack in the ECAL barrel in eta and phi direction More...
 
void doClusteringWorker (const reco::PFRecHitCollection &rechits)
 perform clustering More...
 
void findSeeds (const reco::PFRecHitCollection &rechits)
 look for seeds More...
 
void paint (unsigned rhi, unsigned color=1)
 paint a rechit with a color. More...
 

Private Attributes

bool cleanRBXandHPDs_
 option to clean HCAL RBX's and HPD's More...
 
std::vector< unsigned > color_
 color, for all rechits More...
 
bool debug_
 debugging on/off More...
 
const CaloSubdetectorGeometryeb_geom
 
const CaloSubdetectorGeometryee_geom
 
std::auto_ptr< PositionCalceg_pos_calc
 
std::multimap< double,
unsigned, std::greater< double > > 
eRecHits_
 indices to rechits, sorted by decreasing E (not E_T) More...
 
TFile * file_
 
TH2F * hBNeighbour
 
TH2F * hENeighbour
 
std::set< unsigned > idUsedRecHits_
 ids of rechits used in seed search More...
 
std::vector< bool > mask_
 
std::vector< double > minS4S1Barrel_
 
std::vector< double > minS4S1Endcap_
 
double minS6S2DoubleSpikeBarrel_
 
double minS6S2DoubleSpikeEndcap_
 
int nNeighbours_
 number of neighbours More...
 
double param_W0_
 
std::auto_ptr< std::vector
< reco::PFCluster > > 
pfClusters_
 all clusters More...
 
std::auto_ptr< std::vector
< reco::PFRecHit > > 
pfRecHitsCleaned_
 particle flow rechits cleaned More...
 
int posCalcNCrystal_
 number of crystals for position calculation More...
 
double posCalcP1_
 
const CaloSubdetectorGeometrypreshower_geom
 
PFRecHitHandle rechitsHandle_
 
std::vector< unsigned > seeds_
 vector of indices for seeds. More...
 
std::vector< SeedStateseedStates_
 seed state, for all rechits More...
 
double showerSigma_
 sigma of shower (cm) More...
 
std::auto_ptr
< SortedPFRecHitCollection
sortedRecHits_
 
double threshBarrel_
 barrel threshold More...
 
double threshCleanBarrel_
 Barrel cleaning threshold and S4/S1 smallest fractiom. More...
 
double threshCleanEndcap_
 Endcap cleaning threshold and S4/S1 smallest fractiom. More...
 
double threshDoubleSpikeBarrel_
 Barrel double-spike cleaning. More...
 
double threshDoubleSpikeEndcap_
 Endcap double-spike cleaning. More...
 
double threshEndcap_
 endcap threshold More...
 
double threshPtBarrel_
 
double threshPtEndcap_
 
double threshPtSeedBarrel_
 
double threshPtSeedEndcap_
 
double threshSeedBarrel_
 barrel seed threshold More...
 
double threshSeedEndcap_
 endcap seed threshold More...
 
std::vector< std::vector
< unsigned > > 
topoClusters_
 sets of cells having one common side, and energy over threshold More...
 
bool useCornerCells_
 option to use cells with a common corner to build topo-clusters More...
 
std::vector< bool > usedInTopo_
 used in topo cluster? for all rechits More...
 
PositionCalcType which_pos_calc_
 parameter for position calculation More...
 

Static Private Attributes

static unsigned prodNum_ = 1
 product number More...
 

Friends

std::ostream & operator<< (std::ostream &out, const PFClusterAlgo &algo)
 

Detailed Description

Definition at line 53 of file PFClusterAlgo.h.

Member Typedef Documentation

typedef std::multimap<double, unsigned >::iterator PFClusterAlgo::EH

Definition at line 256 of file PFClusterAlgo.h.

typedef std::map<unsigned, unsigned >::const_iterator PFClusterAlgo::IDH

Definition at line 255 of file PFClusterAlgo.h.

Definition at line 59 of file PFClusterAlgo.h.

Definition at line 79 of file PFClusterAlgo.h.

Definition at line 57 of file PFClusterAlgo.h.

Definition at line 58 of file PFClusterAlgo.h.

Member Enumeration Documentation

Enumerator
NONE 
SEED 
SPECIAL 

Definition at line 207 of file PFClusterAlgo.h.

Returns
threshold, seed threshold, (gaussian width, p1 ??) for a given zone (endcap, barrel, VFCAL ??)
Enumerator
THRESH 
SEED_THRESH 
PT_THRESH 
SEED_PT_THRESH 
CLEAN_THRESH 
CLEAN_S4S1 
DOUBLESPIKE_THRESH 
DOUBLESPIKE_S6S2 

Definition at line 230 of file PFClusterAlgo.h.

Enumerator
EGPositionCalc 
EGPositionFormula 
PFPositionCalc 
kNotDefined 

Definition at line 61 of file PFClusterAlgo.h.

Enumerator
UNKNOWN 
NO 
YES 
CLEAN 

Definition at line 245 of file PFClusterAlgo.h.

Constructor & Destructor Documentation

PFClusterAlgo::PFClusterAlgo ( )

constructor

Definition at line 23 of file PFClusterAlgo.cc.

References file_, hBNeighbour, and hENeighbour.

23  :
24  pfClusters_( new vector<reco::PFCluster> ),
25  pfRecHitsCleaned_( new vector<reco::PFRecHit> ),
26  threshBarrel_(0.),
27  threshPtBarrel_(0.),
28  threshSeedBarrel_(0.2),
30  threshEndcap_(0.),
31  threshPtEndcap_(0.),
32  threshSeedEndcap_(0.6),
34  threshCleanBarrel_(1E5),
35  minS4S1Barrel_(0.),
38  threshCleanEndcap_(1E5),
39  minS4S1Endcap_(0.),
42  nNeighbours_(4),
43  posCalcNCrystal_(-1),
45  posCalcP1_(-1),
46  showerSigma_(5),
47  useCornerCells_(false),
48  cleanRBXandHPDs_(false),
49  debug_(false)
50 {
51  file_ = 0;
52  hBNeighbour = 0;
53  hENeighbour = 0;
54 }
int posCalcNCrystal_
number of crystals for position calculation
double threshEndcap_
endcap threshold
double threshDoubleSpikeBarrel_
Barrel double-spike cleaning.
double threshDoubleSpikeEndcap_
Endcap double-spike cleaning.
double threshPtEndcap_
TH2F * hENeighbour
std::vector< double > minS4S1Barrel_
double threshPtSeedEndcap_
double threshPtSeedBarrel_
std::auto_ptr< std::vector< reco::PFRecHit > > pfRecHitsCleaned_
particle flow rechits cleaned
double threshSeedBarrel_
barrel seed threshold
TH2F * hBNeighbour
bool cleanRBXandHPDs_
option to clean HCAL RBX&#39;s and HPD&#39;s
std::vector< double > minS4S1Endcap_
double threshCleanBarrel_
Barrel cleaning threshold and S4/S1 smallest fractiom.
bool debug_
debugging on/off
int nNeighbours_
number of neighbours
double minS6S2DoubleSpikeEndcap_
PositionCalcType which_pos_calc_
parameter for position calculation
double showerSigma_
sigma of shower (cm)
double threshBarrel_
barrel threshold
bool useCornerCells_
option to use cells with a common corner to build topo-clusters
double threshSeedEndcap_
endcap seed threshold
std::auto_ptr< std::vector< reco::PFCluster > > pfClusters_
all clusters
double threshPtBarrel_
double minS6S2DoubleSpikeBarrel_
double threshCleanEndcap_
Endcap cleaning threshold and S4/S1 smallest fractiom.
virtual PFClusterAlgo::~PFClusterAlgo ( )
inlinevirtual

destructor

Definition at line 70 of file PFClusterAlgo.h.

70 {;}

Member Function Documentation

void PFClusterAlgo::buildPFClusters ( const std::vector< unsigned > &  cluster,
const reco::PFRecHitCollection rechits 
)
private

build PFClusters from a topocluster

Definition at line 1021 of file PFClusterAlgo.cc.

References reco::PFCluster::addRecHitFraction(), calculateClusterPosition(), dtNoiseDBValidation_cfg::cerr, gather_cfg::cout, createRecHitRef(), debug_, delta, diffTreeTool::diff, alignCSCRings::e, relval_parameters_module::energy, create_public_lumi_plots::exp, cropTnPTrees::frac, i, isSeed(), reco::PFRecHit::layer(), max(), paint(), pfClusters_, posCalcNCrystal(), posCalcNCrystal_, reco::PFRecHit::position(), position, rechit(), seedStates_, showerSigma_, SPECIAL, tmp, findQualityFiles::v, and YES.

Referenced by doClusteringWorker().

1023 {
1024 
1025 
1026  // bool debug = false;
1027 
1028 
1029  // several rechits may be seeds. initialize PFClusters on these seeds.
1030 
1031  vector<reco::PFCluster> curpfclusters;
1032  vector<reco::PFCluster> curpfclusterswodepthcor;
1033  vector< unsigned > seedsintopocluster;
1034 
1035 
1036  for(unsigned i=0; i<topocluster.size(); i++ ) {
1037 
1038  unsigned rhi = topocluster[i];
1039 
1040  if( seedStates_[rhi] == YES ) {
1041 
1042  reco::PFCluster cluster;
1043  reco::PFCluster clusterwodepthcor;
1044 
1045  double fraction = 1.0;
1046 
1047  reco::PFRecHitRef recHitRef = createRecHitRef( rechits, rhi );
1048 
1049  cluster.addRecHitFraction( reco::PFRecHitFraction( recHitRef,
1050  fraction ) );
1051 
1052  // cluster.addRecHit( rhi, fraction );
1053 
1054  calculateClusterPosition( cluster,
1055  clusterwodepthcor,
1056  true );
1057 
1058 
1059 // cout<<"PFClusterAlgo: 2"<<endl;
1060  curpfclusters.push_back( cluster );
1061  curpfclusterswodepthcor.push_back( clusterwodepthcor );
1062 #ifdef PFLOW_DEBUG
1063  if(debug_) {
1064  cout << "PFClusterAlgo::buildPFClusters: seed "
1065  << rechit( rhi, rechits) <<endl;
1066  cout << "PFClusterAlgo::buildPFClusters: pfcluster initialized : "
1067  << cluster <<endl;
1068  }
1069 #endif
1070 
1071  // keep track of the seed of each topocluster
1072  seedsintopocluster.push_back( rhi );
1073  }
1074  }
1075 
1076  // if only one seed in the topocluster, use all crystals
1077  // in the position calculation (posCalcNCrystal = -1)
1078  // otherwise, use the user specified value
1079  int posCalcNCrystal = seedsintopocluster.size()>1 ? posCalcNCrystal_:-1;
1080  double ns2 = std::max(1.,(double)(seedsintopocluster.size())-1.);
1081  ns2 *= ns2;
1082 
1083  // Find iteratively the energy and position
1084  // of each pfcluster in the topological cluster
1085  unsigned iter = 0;
1086  unsigned niter = 50;
1087  double diff = ns2;
1088 
1089  // if(debug_) niter=2;
1090  vector<double> ener;
1091  vector<double> dist;
1092  vector<double> frac;
1093  vector<math::XYZVector> tmp;
1094 
1095  while ( iter++ < niter && diff > 1E-8*ns2 ) {
1096 
1097  // Store previous iteration's result and reset pfclusters
1098  ener.clear();
1099  tmp.clear();
1100 
1101  for ( unsigned ic=0; ic<curpfclusters.size(); ic++ ) {
1102  ener.push_back( curpfclusters[ic].energy() );
1103 
1105  v = curpfclusters[ic].position();
1106 
1107  tmp.push_back( v );
1108 
1109 #ifdef PFLOW_DEBUG
1110  if(debug_) {
1111  cout<<"saving photon pos "<<ic<<" "<<curpfclusters[ic]<<endl;
1112  cout<<tmp[ic].X()<<" "<<tmp[ic].Y()<<" "<<tmp[ic].Z()<<endl;
1113  }
1114 #endif
1115 
1116  curpfclusters[ic].reset();
1117  }
1118 
1119  // Loop over topocluster cells
1120  for( unsigned irh=0; irh<topocluster.size(); irh++ ) {
1121  unsigned rhindex = topocluster[irh];
1122 
1123  const reco::PFRecHit& rh = rechit( rhindex, rechits);
1124 
1125  // int layer = rh.layer();
1126 
1127  dist.clear();
1128  frac.clear();
1129  double fractot = 0.;
1130 
1131  bool isaseed = isSeed(rhindex);
1132 
1133  math::XYZVector cposxyzcell;
1134  cposxyzcell = rh.position();
1135 
1136 #ifdef PFLOW_DEBUG
1137  if(debug_) {
1138  cout<<rh<<endl;
1139  cout<<"start loop on curpfclusters"<<endl;
1140  }
1141 #endif
1142 
1143  // Loop over pfclusters
1144  for ( unsigned ic=0; ic<tmp.size(); ic++) {
1145 
1146 #ifdef PFLOW_DEBUG
1147  if(debug_) cout<<"pfcluster "<<ic<<endl;
1148 #endif
1149 
1150  double frc=0.;
1151  bool seedexclusion=true;
1152 
1153  // convert cluster coordinates to xyz
1154  //math::XYZVector cposxyzclust( tmp[ic].X(), tmp[ic].Y(), tmp[ic].Z() );
1155  // cluster position used to compute distance with cell
1156  math::XYZVector cposxyzclust;
1157  cposxyzclust = curpfclusterswodepthcor[ic].position();
1158 
1159 #ifdef PFLOW_DEBUG
1160  if(debug_) {
1161 
1162  cout<<"CLUSTER "<<cposxyzclust.X()<<","
1163  <<cposxyzclust.Y()<<","
1164  <<cposxyzclust.Z()<<"\t\t"
1165  <<"CELL "<<cposxyzcell.X()<<","
1166  <<cposxyzcell.Y()<<","
1167  <<cposxyzcell.Z()<<endl;
1168  }
1169 #endif
1170 
1171  // Compute the distance between the current cell
1172  // and the current PF cluster, normalized to a
1173  // number of "sigma"
1174  math::XYZVector deltav = cposxyzclust;
1175  deltav -= cposxyzcell;
1176  double d = deltav.R() / showerSigma_;
1177 
1178  // if distance cell-cluster is too large, it means that
1179  // we're right on the junction between 2 subdetectors (HCAL/VFCAL...)
1180  // in this case, distance is calculated in the xy plane
1181  // could also be a large supercluster...
1182 #ifdef PFLOW_DEBUG
1183  if( d > 10. && debug_ ) {
1184  paint(rhindex, SPECIAL);
1185  cout<<"PFClusterAlgo Warning: distance too large"<<d<<endl;
1186  }
1187 #endif
1188  dist.push_back( d );
1189 
1190  // the current cell is the seed from the current photon.
1191  if( rhindex == seedsintopocluster[ic] && seedexclusion ) {
1192  frc = 1.;
1193 #ifdef PFLOW_DEBUG
1194  if(debug_) cout<<"this cell is a seed for the current photon"<<endl;
1195 #endif
1196  }
1197  else if( isaseed && seedexclusion ) {
1198  frc = 0.;
1199 #ifdef PFLOW_DEBUG
1200  if(debug_) cout<<"this cell is a seed for another photon"<<endl;
1201 #endif
1202  }
1203  else {
1204  // Compute the fractions of the cell energy to be assigned to
1205  // each curpfclusters in the cluster.
1206  frc = ener[ic] * exp ( - dist[ic]*dist[ic] / 2. );
1207 
1208 #ifdef PFLOW_DEBUG
1209  if(debug_) {
1210  cout<<"dist["<<ic<<"] "<<dist[ic]
1211  // <<", sigma="<<sigma
1212  <<", frc="<<frc<<endl;
1213  }
1214 #endif
1215 
1216  }
1217  fractot += frc;
1218  frac.push_back(frc);
1219  }
1220 
1221  // Add the relevant fraction of the cell to the curpfclusters
1222 #ifdef PFLOW_DEBUG
1223  if(debug_) cout<<"start add cell"<<endl;
1224 #endif
1225  for ( unsigned ic=0; ic<tmp.size(); ++ic ) {
1226 #ifdef PFLOW_DEBUG
1227  if(debug_)
1228  cout<<" frac["<<ic<<"] "<<frac[ic]<<" "<<fractot<<" "<<rh<<endl;
1229 #endif
1230 
1231  if( fractot )
1232  frac[ic] /= fractot;
1233  else {
1234 #ifdef PFLOW_DEBUG
1235  if( debug_ ) {
1236  int layer = rh.layer();
1237  cerr<<"fractot = 0 ! "<<layer<<endl;
1238 
1239  for( unsigned trh=0; trh<topocluster.size(); trh++ ) {
1240  unsigned tindex = topocluster[trh];
1241  const reco::PFRecHit& rh = rechit( tindex, rechits);
1242  cout<<rh<<endl;
1243  }
1244 
1245  // assert(0)
1246  }
1247 #endif
1248 
1249  continue;
1250  }
1251 
1252  // if the fraction has been set to 0, the cell
1253  // is now added to the cluster - careful ! (PJ, 19/07/08)
1254  // BUT KEEP ONLY CLOSE CELLS OTHERWISE MEMORY JUST EXPLOSES
1255  // (PJ, 15/09/08 <- similar to what existed before the
1256  // previous bug fix, but keeps the close seeds inside,
1257  // even if their fraction was set to zero.)
1258  // Also add a protection to keep the seed in the cluster
1259  // when the latter gets far from the former. These cases
1260  // (about 1% of the clusters) need to be studied, as
1261  // they create fake photons, in general.
1262  // (PJ, 16/09/08)
1263  if ( dist[ic] < 10. || frac[ic] > 0.99999 ) {
1264  // if ( dist[ic] > 6. ) cout << "Warning : PCluster is getting very far from its seeding cell" << endl;
1265  reco::PFRecHitRef recHitRef = createRecHitRef( rechits, rhindex );
1266  reco::PFRecHitFraction rhf( recHitRef,frac[ic] );
1267  curpfclusters[ic].addRecHitFraction( rhf );
1268  }
1269  }
1270  // if(debug_) cout<<" end add cell"<<endl;
1271  }
1272 
1273  // Determine the new cluster position and check
1274  // the distance with the previous iteration
1275  diff = 0.;
1276  for ( unsigned ic=0; ic<tmp.size(); ++ic ) {
1277 
1278  calculateClusterPosition( curpfclusters[ic], curpfclusterswodepthcor[ic],
1279  true, posCalcNCrystal );
1280 #ifdef PFLOW_DEBUG
1281  if(debug_) cout<<"new iter "<<ic<<endl;
1282  if(debug_) cout<<curpfclusters[ic]<<endl;
1283 #endif
1284 
1285  double delta = ROOT::Math::VectorUtil::DeltaR(curpfclusters[ic].position(),tmp[ic]);
1286  if ( delta > diff ) diff = delta;
1287  }
1288  }
1289 
1290  // Issue a warning message if the number of iterations
1291  // exceeds 50
1292 #ifdef PFLOW_DEBUG
1293  if ( iter >= 50 && debug_ )
1294  cout << "PFClusterAlgo Warning: "
1295  << "more than "<<niter<<" iterations in pfcluster finding: "
1296  << setprecision(10) << diff << endl;
1297 #endif
1298 
1299  // There we go
1300  // add all clusters to the list of pfClusters.
1301  // But first clean the clusters so that rechits with negligible fraction
1302  // are removed.
1303  for(unsigned ic=0; ic<curpfclusters.size(); ic++) {
1304 
1305  //copy full list of RecHitFractions
1306  std::vector< reco::PFRecHitFraction > rhfracs = curpfclusters[ic].recHitFractions();
1307  //reset cluster
1308  curpfclusters[ic].reset();
1309 
1310  //loop over full list of rechit fractions and add the back to the cluster only
1311  //if the fraction is above some reasonable threshold
1312  for (const auto &rhf : rhfracs) {
1313  if (rhf.fraction()>1e-7) {
1314  curpfclusters[ic].addRecHitFraction(rhf);
1315  }
1316  }
1317 
1318  calculateClusterPosition(curpfclusters[ic], curpfclusterswodepthcor[ic],
1319  true, posCalcNCrystal);
1320 
1321  pfClusters_->push_back(curpfclusters[ic]);
1322  }
1323 }
dbl * delta
Definition: mlp_gen.cc:36
int i
Definition: DBlmapReader.cc:9
int posCalcNCrystal_
number of crystals for position calculation
Particle flow cluster, see clustering algorithm in PFClusterAlgo.
Definition: PFCluster.h:43
bool isSeed(unsigned rhi) const
const reco::PFRecHit & rechit(unsigned i, const reco::PFRecHitCollection &rechits)
Fraction of a PFRecHit (rechits can be shared between several PFCluster&#39;s)
std::vector< SeedState > seedStates_
seed state, for all rechits
static int position[TOTALCHAMBERS][3]
Definition: ReadPGInfo.cc:509
PFLayer::Layer layer() const
rechit layer
Definition: PFRecHit.h:108
int posCalcNCrystal() const
get number of crystals for position calculation (-1 all,5, or 9)
Particle flow rechit (rechit + geometry and topology information). See clustering algorithm in PFClus...
Definition: PFRecHit.h:31
const T & max(const T &a, const T &b)
const math::XYZPoint & position() const
is seed ?
Definition: PFRecHit.h:141
bool debug_
debugging on/off
reco::PFRecHitRef createRecHitRef(const reco::PFRecHitCollection &rechits, unsigned rhi)
XYZVectorD XYZVector
spatial vector with cartesian internal representation
Definition: Vector3D.h:30
double showerSigma_
sigma of shower (cm)
void paint(unsigned rhi, unsigned color=1)
paint a rechit with a color.
void addRecHitFraction(const reco::PFRecHitFraction &frac)
add a given fraction of the rechit
Definition: PFCluster.cc:40
std::vector< std::vector< double > > tmp
Definition: MVATrainer.cc:100
void calculateClusterPosition(reco::PFCluster &cluster, reco::PFCluster &clusterwodepthcor, bool depcor=true, int posCalcNCrystal=0)
calculate position of a cluster
tuple cout
Definition: gather_cfg.py:121
std::auto_ptr< std::vector< reco::PFCluster > > pfClusters_
all clusters
void PFClusterAlgo::buildTopoCluster ( std::vector< unsigned > &  cluster,
unsigned  rhi,
const reco::PFRecHitCollection rechits 
)
private

build a topocluster (recursive)

Definition at line 945 of file PFClusterAlgo.cc.

References funct::abs(), gather_cfg::cout, debug_, alignCSCRings::e, reco::PFRecHit::energy(), PFLayer::HCAL_BARREL2, i, reco::PFRecHit::layer(), masked(), reco::PFRecHit::neighbours4(), reco::PFRecHit::neighbours8(), parameter(), reco::PFRecHit::positionREP(), reco::PFRecHit::pt2(), PT_THRESH, rechit(), GOODCOLL_filter_cfg::thresh, THRESH, useCornerCells_, and usedInTopo_.

Referenced by buildTopoClusters().

947  {
948 
949 
950 #ifdef PFLOW_DEBUG
951  if(debug_)
952  cout<<"PFClusterAlgo::buildTopoCluster in"<<endl;
953 #endif
954 
955  const reco::PFRecHit& rh = rechit( rhi, rechits);
956 
957  double e = rh.energy();
958  int layer = rh.layer();
959 
960 
961  int iring = 0;
962  if (layer==PFLayer::HCAL_BARREL2 && abs(rh.positionREP().Eta())>0.34) iring= 1;
963 
964  double thresh = parameter( THRESH,
965  static_cast<PFLayer::Layer>(layer), 0, iring );
966  double ptThresh = parameter( PT_THRESH,
967  static_cast<PFLayer::Layer>(layer), 0, iring );
968 
969 
970  if( e < thresh || (ptThresh > 0. && rh.pt2() < ptThresh*ptThresh) ) {
971 #ifdef PFLOW_DEBUG
972  if(debug_)
973  cout<<"return : "<<e<<"<"<<thresh<<endl;
974 #endif
975  return;
976  }
977 
978  // add hit to cluster
979 
980  cluster.push_back( rhi );
981  // idUsedRecHits_.insert( rh.detId() );
982 
983  usedInTopo_[ rhi ] = true;
984 
985  // cout<<" hit ptr "<<hit<<endl;
986 
987  // get neighbours, either with one side in common,
988  // or with one corner in common (if useCornerCells_)
989  const std::vector< unsigned >& nbs
990  = useCornerCells_ ? rh.neighbours8() : rh.neighbours4();
991 
992  for(unsigned i=0; i<nbs.size(); i++) {
993 
994 // const reco::PFRecHit& neighbour = rechit( nbs[i], rechits );
995 
996 // set<unsigned>::iterator used
997 // = idUsedRecHits_.find( neighbour.detId() );
998 // if(used != idUsedRecHits_.end() ) continue;
999 
1000  // already used
1001  if( usedInTopo_[ nbs[i] ] ) {
1002 #ifdef PFLOW_DEBUG
1003  if(debug_)
1004  cout<<rhi<<" used"<<endl;
1005 #endif
1006  continue;
1007  }
1008 
1009  if( !masked(nbs[i]) ) continue;
1010  buildTopoCluster( cluster, nbs[i], rechits );
1011  }
1012 #ifdef PFLOW_DEBUG
1013  if(debug_)
1014  cout<<"PFClusterAlgo::buildTopoCluster out"<<endl;
1015 #endif
1016 
1017 }
int i
Definition: DBlmapReader.cc:9
double pt2() const
rechit momentum transverse to the beam, squared.
Definition: PFRecHit.h:123
const reco::PFRecHit & rechit(unsigned i, const reco::PFRecHitCollection &rechits)
std::vector< bool > usedInTopo_
used in topo cluster? for all rechits
PFLayer::Layer layer() const
rechit layer
Definition: PFRecHit.h:108
bool masked(unsigned rhi) const
Particle flow rechit (rechit + geometry and topology information). See clustering algorithm in PFClus...
Definition: PFRecHit.h:31
const std::vector< unsigned > & neighbours8() const
Definition: PFRecHit.h:160
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
bool debug_
debugging on/off
void buildTopoCluster(std::vector< unsigned > &cluster, unsigned rhi, const reco::PFRecHitCollection &rechits)
build a topocluster (recursive)
double parameter(Parameter paramtype, PFLayer::Layer layer, unsigned iCoeff=0, int iring0=0) const
double energy() const
rechit energy
Definition: PFRecHit.h:111
bool useCornerCells_
option to use cells with a common corner to build topo-clusters
tuple cout
Definition: gather_cfg.py:121
const std::vector< unsigned > & neighbours4() const
Definition: PFRecHit.h:157
const REPPoint & positionREP() const
rechit cell centre rho, eta, phi. call calculatePositionREP before !
Definition: PFRecHit.cc:121
void PFClusterAlgo::buildTopoClusters ( const reco::PFRecHitCollection rechits)
private

build topoclusters around seeds

Definition at line 901 of file PFClusterAlgo.cc.

References buildTopoCluster(), gather_cfg::cout, debug_, masked(), seeds_, topoClusters_, and usedInTopo_.

Referenced by doClusteringWorker().

901  {
902 
903  topoClusters_.clear();
904 
905 #ifdef PFLOW_DEBUG
906  if(debug_)
907  cout<<"PFClusterAlgo::buildTopoClusters start"<<endl;
908 #endif
909 
910  for(unsigned is = 0; is<seeds_.size(); is++) {
911 
912  unsigned rhi = seeds_[is];
913 
914  if( !masked(rhi) ) continue;
915  // rechit was masked to be processed
916 
917  // already used in a topological cluster
918  if( usedInTopo_[rhi] ) {
919 #ifdef PFLOW_DEBUG
920  if(debug_)
921  cout<<rhi<<" used"<<endl;
922 #endif
923  continue;
924  }
925 
926  vector< unsigned > topocluster;
927  buildTopoCluster( topocluster, rhi, rechits );
928 
929  if(topocluster.empty() ) continue;
930 
931  topoClusters_.push_back( topocluster );
932 
933  }
934 
935 #ifdef PFLOW_DEBUG
936  if(debug_)
937  cout<<"PFClusterAlgo::buildTopoClusters done"<<endl;
938 #endif
939 
940  return;
941 }
std::vector< std::vector< unsigned > > topoClusters_
sets of cells having one common side, and energy over threshold
std::vector< bool > usedInTopo_
used in topo cluster? for all rechits
bool masked(unsigned rhi) const
std::vector< unsigned > seeds_
vector of indices for seeds.
bool debug_
debugging on/off
void buildTopoCluster(std::vector< unsigned > &cluster, unsigned rhi, const reco::PFRecHitCollection &rechits)
build a topocluster (recursive)
tuple cout
Definition: gather_cfg.py:121
void PFClusterAlgo::calculateClusterPosition ( reco::PFCluster cluster,
reco::PFCluster clusterwodepthcor,
bool  depcor = true,
int  posCalcNCrystal = 0 
)
private

calculate position of a cluster

Definition at line 1328 of file PFClusterAlgo.cc.

References funct::abs(), dtNoiseDBValidation_cfg::cerr, gather_cfg::cout, reco::PFCluster::depthCorA_, reco::PFCluster::depthCorAp_, reco::PFCluster::depthCorB_, reco::PFCluster::depthCorBp_, reco::PFCluster::depthCorMode_, reco::PFRecHit::detId(), alignCSCRings::e, eb_geom, PFLayer::ECAL_BARREL, PFLayer::ECAL_ENDCAP, ee_geom, eg_pos_calc, EGPositionCalc, EGPositionFormula, reco::PFCluster::energy(), reco::PFRecHit::energy(), reco::CaloCluster::energy_, edm::hlt::Exception, reco::PFRecHit::getAxisXYZ(), PFLayer::HCAL_BARREL1, PFLayer::HCAL_BARREL2, PFLayer::HCAL_ENDCAP, PFLayer::HF_EM, PFLayer::HF_HAD, reco::PFRecHit::isNeighbour4(), reco::PFRecHit::isNeighbour8(), isSeed(), reco::PFCluster::layer(), reco::PFRecHit::layer(), LayerTriplets::layers(), create_public_lumi_plots::log, max(), PFLayer::NONE, SVfit_namespace::normalize(), NULL, p1, param_W0_, PFPositionCalc, posCalcNCrystal_, posCalcP1_, reco::PFRecHit::position(), reco::CaloCluster::position_, reco::PFCluster::posrep_, preshower_geom, PFLayer::PS1, PFLayer::PS2, reco::PFCluster::rechits_, reco::PFCluster::setLayer(), reco::CaloCluster::setSeed(), sortedRecHits_, mathSSE::sqrt(), threshBarrel_, threshEndcap_, which_pos_calc_, x, detailsBasic3DVector::y, and detailsBasic3DVector::z.

Referenced by buildPFClusters().

1331  {
1332 
1333  if( posCalcNCrystal_ != -1 &&
1334  posCalcNCrystal_ != 5 &&
1335  posCalcNCrystal_ != 9 ) {
1336  throw "PFCluster::calculatePosition : posCalcNCrystal_ must be -1, 5, or 9.";
1337  }
1338 
1339  std::vector< std::pair< DetId, float > > hits_and_fracts;
1340 
1342 
1343  cluster.position_.SetXYZ(0,0,0);
1344 
1345  cluster.energy_ = 0;
1346 
1347  double normalize = 0;
1348 
1349  // calculate total energy, average layer, and look for seed ---------- //
1350 
1351 
1352  // double layer = 0;
1353  map <PFLayer::Layer, double> layers;
1354 
1355  unsigned seedIndex = 0;
1356  bool seedIndexFound = false;
1357 
1358  //Colin: the following can be simplified!
1359 
1360  // loop on rechit fractions
1361  for (unsigned ic=0; ic<cluster.rechits_.size(); ic++ ) {
1362 
1363  unsigned rhi = cluster.rechits_[ic].recHitRef().index();
1364  // const reco::PFRecHit& rh = rechit( rhi, rechits );
1365 
1366  const reco::PFRecHit& rh = *(cluster.rechits_[ic].recHitRef());
1367  double fraction = cluster.rechits_[ic].fraction();
1368 
1369  // Find the seed of this sub-cluster (excluding other seeds found in the topological
1370  // cluster, the energy fraction of which were set to 0 fpr the position determination.
1371  if( isSeed(rhi) && fraction > 1e-9 ) {
1372  seedIndex = rhi;
1373  seedIndexFound = true;
1374  cluster.setSeed(DetId(cluster.rechits_[ic].recHitRef()->detId()));
1375  }
1376 
1377  double recHitEnergy = rh.energy() * fraction;
1378 
1379  // is nan ?
1380  if( recHitEnergy!=recHitEnergy ) {
1381  ostringstream ostr;
1382  edm::LogError("PFClusterAlgo")<<"rechit "<<rh.detId()<<" has a NaN energy... The input of the particle flow clustering seems to be corrupted.";
1383  }
1384 
1385  cluster.energy_ += recHitEnergy;
1386 
1387  // sum energy in each layer
1388  PFLayer::Layer layer = rh.layer();
1389 
1390  map <PFLayer::Layer, double>:: iterator it = layers.find(layer);
1391 
1392  if (it != layers.end()) {
1393  it->second += recHitEnergy;
1394  } else {
1395 
1396  layers.insert(make_pair(layer, recHitEnergy));
1397 
1398  }
1399 
1400  }
1401 
1402  assert(seedIndexFound);
1403 
1404  // loop over pairs to find layer with max energy
1405  double Emax = 0.;
1406  PFLayer::Layer layer = PFLayer::NONE;
1407  for (map<PFLayer::Layer, double>::iterator it = layers.begin();
1408  it != layers.end(); ++it) {
1409  double e = it->second;
1410  if(e > Emax){
1411  Emax = e;
1412  layer = it->first;
1413  }
1414  }
1415 
1416 
1417  //setlayer here
1418  cluster.setLayer( layer ); // take layer with max energy
1419 
1420  // layer /= cluster.energy_;
1421  // cluster.layer_ = lrintf(layer); // nearest integer
1422 
1423 
1424 
1425  double p1 = posCalcP1_;
1426  if( p1 < 0 ) {
1427  // automatic (and hopefully best !) determination of the parameter
1428  // for position determination.
1429 
1430  // Remove the ad-hoc determination of p1, and set it to the
1431  // seed threshold.
1432  switch(cluster.layer() ) {
1433  case PFLayer::ECAL_BARREL:
1434  case PFLayer::HCAL_BARREL1:
1435  case PFLayer::HCAL_BARREL2:
1436  // case PFLayer::HCAL_HO:
1437  p1 = threshBarrel_;
1438  break;
1439  case PFLayer::ECAL_ENDCAP:
1440  case PFLayer::HCAL_ENDCAP:
1441  case PFLayer::HF_EM:
1442  case PFLayer::HF_HAD:
1443  case PFLayer::PS1:
1444  case PFLayer::PS2:
1445  p1 = threshEndcap_;
1446  break;
1447 
1448  /*
1449  switch(cluster.layer() ) {
1450  case PFLayer::ECAL_BARREL:
1451  p1 = 0.004 + 0.022*cluster.energy_; // 27 feb 2006
1452  break;
1453  case PFLayer::ECAL_ENDCAP:
1454  p1 = 0.014 + 0.009*cluster.energy_; // 27 feb 2006
1455  break;
1456  case PFLayer::HCAL_BARREL1:
1457  case PFLayer::HCAL_BARREL2:
1458  case PFLayer::HCAL_ENDCAP:
1459  case PFLayer::HCAL_HF:
1460  p1 = 5.41215e-01 * log( cluster.energy_ / 1.29803e+01 );
1461  if(p1<0.01) p1 = 0.01;
1462  break;
1463  */
1464 
1465  default:
1466  cerr<<"Clusters weight_p1 -1 not yet allowed for layer "<<layer
1467  <<". Chose a better value in the opt file"<<endl;
1468  assert(0);
1469  break;
1470  }
1471  }
1472  else if( p1< 1e-9 ) { // will divide by p1 later on
1473  p1 = 1e-9;
1474  }
1475  // calculate uncorrected cluster position --------------------------------
1476 
1477  reco::PFCluster::REPPoint clusterpos; // uncorrected cluster pos
1478  math::XYZPoint clusterposxyz; // idem, xyz coord
1479  math::XYZPoint firstrechitposxyz; // pos of the rechit with highest E
1480 
1481  double maxe = -9999;
1482  double x = 0;
1483  double y = 0;
1484  double z = 0;
1485 
1486  for (unsigned ic=0; ic<cluster.rechits_.size(); ic++ ) {
1487 
1488  unsigned rhi = cluster.rechits_[ic].recHitRef().index();
1489 // const reco::PFRecHit& rh = rechit( rhi, rechits );
1490 
1491  const reco::PFRecHit& rh = *(cluster.rechits_[ic].recHitRef());
1492 
1493  hits_and_fracts.push_back(std::make_pair(DetId(rh.detId()),
1494  cluster.rechits_[ic].fraction()));
1495  if(rhi != seedIndex) { // not the seed
1496  if( posCalcNCrystal == 5 ) { // pos calculated from the 5 neighbours only
1497  if(!rh.isNeighbour4(seedIndex) ) {
1498  continue;
1499  }
1500  }
1501  if( posCalcNCrystal == 9 ) { // pos calculated from the 9 neighbours only
1502  if(!rh.isNeighbour8(seedIndex) ) {
1503  continue;
1504  }
1505  }
1506  }
1507 
1508  double fraction = hits_and_fracts.back().second;
1509  double recHitEnergy = rh.energy() * fraction;
1510  double norm = fraction < 1E-9 ? 0. : max(0., log(recHitEnergy/p1));
1511 
1512  const math::XYZPoint& rechitposxyz = rh.position();
1513 
1514  if( recHitEnergy > maxe ) {
1515  firstrechitposxyz = rechitposxyz;
1516  maxe = recHitEnergy;
1517  }
1518 
1519  x += rechitposxyz.X() * norm;
1520  y += rechitposxyz.Y() * norm;
1521  z += rechitposxyz.Z() * norm;
1522 
1523  // clusterposxyz += rechitposxyz * norm;
1524  normalize += norm;
1525  }
1526 
1527  // normalize uncorrected position
1528  // assert(normalize);
1529  if( normalize < 1e-9 ) {
1530  // cerr<<"--------------------"<<endl;
1531  // cerr<<(*this)<<endl;
1532  cout << "Watch out : cluster too far from its seeding cell, set to 0,0,0" << endl;
1533  clusterposxyz.SetXYZ(0,0,0);
1534  clusterpos.SetCoordinates(0,0,0);
1535  return;
1536  }
1537  else {
1538  x /= normalize;
1539  y /= normalize;
1540  z /= normalize;
1541 
1542  clusterposxyz.SetCoordinates( x, y, z);
1543 
1544  clusterpos.SetCoordinates( clusterposxyz.Rho(), clusterposxyz.Eta(), clusterposxyz.Phi() );
1545 
1546  }
1547 
1548  cluster.posrep_ = clusterpos;
1549 
1550  cluster.position_ = clusterposxyz;
1551 
1552  clusterwodepthcor = cluster;
1553 
1554 
1555  // correctio of the rechit position,
1556  // according to the depth, only for ECAL
1557 
1558 
1559  if( depcor && // correction requested and ECAL
1560  ( cluster.layer() == PFLayer::ECAL_BARREL ||
1561  cluster.layer() == PFLayer::ECAL_ENDCAP ) ) {
1562  if( which_pos_calc_ == EGPositionCalc ) {
1563  // calculate using EG position calc directly
1564  math::XYZPoint clusterposxyzcor(0,0,0);
1565  switch(cluster.layer()) {
1566  case PFLayer::ECAL_BARREL:
1567  clusterposxyzcor = eg_pos_calc->Calculate_Location(hits_and_fracts,
1568  sortedRecHits_.get(),
1569  eb_geom,
1570  NULL);
1571  break;
1572  case PFLayer::ECAL_ENDCAP:
1573  clusterposxyzcor = eg_pos_calc->Calculate_Location(hits_and_fracts,
1574  sortedRecHits_.get(),
1575  ee_geom,
1576  preshower_geom);
1577  break;
1578  default:
1579  break;
1580  }
1581  cluster.posrep_.SetCoordinates( clusterposxyzcor.Rho(),
1582  clusterposxyzcor.Eta(),
1583  clusterposxyzcor.Phi() );
1584  cluster.position_ = clusterposxyzcor;
1585  clusterposxyz = clusterposxyzcor;
1586 
1587  } else { // allow position formula or PF formula
1588 
1589  double corra = reco::PFCluster::depthCorA_;
1590  double corrb = reco::PFCluster::depthCorB_;
1591  if( abs(clusterpos.Eta() ) < 2.6 &&
1592  abs(clusterpos.Eta() ) > 1.65 ) {
1593  // if crystals under preshower, correction is not the same
1594  // (shower depth smaller)
1597  }
1598 
1599  double depth = 0;
1600 
1601  switch( reco::PFCluster::depthCorMode_ ) {
1602  case 1: // for e/gamma
1603  depth = corra * ( corrb + log(cluster.energy_) );
1604  break;
1605  case 2: // for hadrons
1606  depth = corra;
1607  break;
1608  default:
1609  throw cms::Exception("InvalidOption")
1610  <<"PFClusterAlgo::calculateClusterPosition : unknown function"
1611  <<" for depth correction! "<<endl;
1612  }
1613 
1614 
1615  // calculate depth vector:
1616  // its mag is depth
1617  // its direction is the cluster direction (uncorrected)
1618 
1619  // double xdepthv = clusterposxyz.X();
1620  // double ydepthv = clusterposxyz.Y();
1621  // double zdepthv = clusterposxyz.Z();
1622  // double mag = sqrt( xdepthv*xdepthv +
1623  // ydepthv*ydepthv +
1624  // zdepthv*zdepthv );
1625 
1626 
1627  // math::XYZPoint depthv(clusterposxyz);
1628  // depthv.SetMag(depth);
1629 
1630 
1631  math::XYZVector depthv( clusterposxyz.X(),
1632  clusterposxyz.Y(),
1633  clusterposxyz.Z() );
1634  depthv /= sqrt(depthv.Mag2() );
1635  depthv *= depth;
1636 
1637 
1638  // now calculate corrected cluster position:
1639  math::XYZPoint clusterposxyzcor;
1640 
1641  maxe = -9999;
1642  x = 0;
1643  y = 0;
1644  z = 0;
1645  cluster.posrep_.SetXYZ(0,0,0);
1646  normalize = 0;
1647 
1648  for (unsigned ic=0; ic<cluster.rechits_.size(); ic++ ) {
1649 
1650  unsigned rhi = cluster.rechits_[ic].recHitRef().index();
1651  // const reco::PFRecHit& rh = rechit( rhi, rechits );
1652 
1653  const reco::PFRecHit& rh = *(cluster.rechits_[ic].recHitRef());
1654 
1655  if(rhi != seedIndex) {
1656  if( posCalcNCrystal == 5 ) {
1657  if(!rh.isNeighbour4(seedIndex) ) {
1658  continue;
1659  }
1660  }
1661  if( posCalcNCrystal == 9 ) {
1662  if(!rh.isNeighbour8(seedIndex) ) {
1663  continue;
1664  }
1665  }
1666  }
1667 
1668  double fraction = cluster.rechits_[ic].fraction();
1669  double recHitEnergy = rh.energy() * fraction;
1670 
1671  const math::XYZPoint& rechitposxyz = rh.position();
1672 
1673  // rechit axis not correct !
1674  math::XYZVector rechitaxis = rh.getAxisXYZ();
1675  // rechitaxis -= math::XYZVector( rechitposxyz.X(), rechitposxyz.Y(), rechitposxyz.Z() );
1676 
1677  math::XYZVector rechitaxisu( rechitaxis );
1678  rechitaxisu /= sqrt( rechitaxis.Mag2() );
1679 
1680  math::XYZVector displacement( rechitaxisu );
1681  // displacement /= sqrt( displacement.Mag2() );
1682  displacement *= rechitaxisu.Dot( depthv );
1683 
1684  math::XYZPoint rechitposxyzcor( rechitposxyz );
1685  rechitposxyzcor += displacement;
1686 
1687  if( recHitEnergy > maxe ) {
1688  firstrechitposxyz = rechitposxyzcor;
1689  maxe = recHitEnergy;
1690  }
1691 
1692  double norm = -1;
1693  double log_efrac = -1.0;
1694  switch(which_pos_calc_) {
1695  case EGPositionCalc: // in case strange things are happening
1696  case EGPositionFormula:
1697  if( recHitEnergy > 0 ) {
1698  log_efrac = std::log(recHitEnergy/cluster.energy());
1699  }
1700  norm = (log_efrac != -1.0)*std::max(0.0, param_W0_ + log_efrac);
1701  break;
1702  case PFPositionCalc:
1703  norm = fraction < 1E-9 ? 0. : max(0., log(recHitEnergy/p1));
1704  break;
1705  default:
1706  throw cms::Exception("InvalidOption")
1707  << "Invalid position calc type chosen: " << which_pos_calc_ << "!";
1708  }
1709 
1710  x += rechitposxyzcor.X() * norm;
1711  y += rechitposxyzcor.Y() * norm;
1712  z += rechitposxyzcor.Z() * norm;
1713 
1714  // clusterposxyzcor += rechitposxyzcor * norm;
1715  normalize += norm;
1716  }
1717  // normalize
1718  if(normalize < 1e-9) {
1719  cerr<<"--------------------"<<endl;
1720  cerr<< cluster <<endl;
1721  assert(0);
1722  } else {
1723  x /= normalize;
1724  y /= normalize;
1725  z /= normalize;
1726 
1727 
1728  clusterposxyzcor.SetCoordinates(x,y,z);
1729  cluster.posrep_.SetCoordinates( clusterposxyzcor.Rho(),
1730  clusterposxyzcor.Eta(),
1731  clusterposxyzcor.Phi() );
1732  cluster.position_ = clusterposxyzcor;
1733  clusterposxyz = clusterposxyzcor;
1734  }
1735  }
1736  }
1737 }
PFLayer::Layer layer() const
cluster layer, see PFLayer.h in this directory
Definition: PFCluster.cc:81
void setLayer(PFLayer::Layer layer)
set layer
Definition: PFCluster.cc:74
math::XYZPoint position_
cluster centroid position
Definition: CaloCluster.h:209
int posCalcNCrystal_
number of crystals for position calculation
double threshEndcap_
endcap threshold
std::vector< LayerSetAndLayers > layers(const SeedingLayerSetsHits &sets)
Definition: LayerTriplets.cc:4
double energy_
cluster energy
Definition: CaloCluster.h:205
const CaloSubdetectorGeometry * eb_geom
unsigned detId() const
rechit detId
Definition: PFRecHit.h:105
std::vector< reco::PFRecHitFraction > rechits_
vector of rechit fractions (transient)
Definition: PFCluster.h:140
ROOT::Math::PositionVector3D< ROOT::Math::CylindricalEta3D< Double32_t > > REPPoint
Definition: PFCluster.h:47
bool isSeed(unsigned rhi) const
#define NULL
Definition: scimark2.h:8
float float float z
void setSeed(const DetId &id)
Definition: CaloCluster.h:117
PFLayer::Layer layer() const
rechit layer
Definition: PFRecHit.h:108
int posCalcNCrystal() const
get number of crystals for position calculation (-1 all,5, or 9)
Particle flow rechit (rechit + geometry and topology information). See clustering algorithm in PFClus...
Definition: PFRecHit.h:31
const T & max(const T &a, const T &b)
const CaloSubdetectorGeometry * preshower_geom
T sqrt(T t)
Definition: SSEVec.h:48
const math::XYZVector & getAxisXYZ() const
rechit cell axis x, y, z
Definition: PFRecHit.h:147
static double depthCorAp_
Definition: PFCluster.h:156
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
reco::Candidate::Vector normalize(const reco::Candidate::Vector &p)
std::auto_ptr< PositionCalc > eg_pos_calc
static double depthCorB_
Definition: PFCluster.h:153
static double depthCorBp_
Definition: PFCluster.h:159
double energy() const
cluster energy
Definition: PFCluster.h:72
Layer
layer definition
Definition: PFLayer.h:31
static int depthCorMode_
Definition: PFCluster.h:147
const math::XYZPoint & position() const
is seed ?
Definition: PFRecHit.h:141
Definition: DetId.h:18
REPPoint posrep_
cluster position: rho, eta, phi (transient)
Definition: PFCluster.h:143
XYZVectorD XYZVector
spatial vector with cartesian internal representation
Definition: Vector3D.h:30
bool isNeighbour4(unsigned id) const
Definition: PFRecHit.cc:248
XYZPointD XYZPoint
point in space with cartesian internal representation
Definition: Point3D.h:12
PositionCalcType which_pos_calc_
parameter for position calculation
double energy() const
rechit energy
Definition: PFRecHit.h:111
double p1[4]
Definition: TauolaWrapper.h:89
double threshBarrel_
barrel threshold
bool isNeighbour8(unsigned id) const
Definition: PFRecHit.cc:257
const CaloSubdetectorGeometry * ee_geom
tuple cout
Definition: gather_cfg.py:121
std::auto_ptr< SortedPFRecHitCollection > sortedRecHits_
static double depthCorA_
Definition: PFCluster.h:150
Definition: DDAxes.h:10
void PFClusterAlgo::cleanRBXAndHPD ( const reco::PFRecHitCollection rechits)
private

Clean HCAL readout box noise and HPD discharge.

Definition at line 309 of file PFClusterAlgo.cc.

References funct::abs(), reco::PFRecHit::detId(), reco::PFRecHit::energy(), eRecHits_, PFLayer::HCAL_BARREL1, PFLayer::HCAL_ENDCAP, HcalDetId::ieta(), recoMuon::in, HcalDetId::iphi(), reco::PFRecHit::layer(), python.multivaluedict::map(), mask_, masked(), reco::PFRecHit::neighbours4(), pfRecHitsCleaned_, reco::PFRecHit::position(), rechit(), mathSSE::sqrt(), and dtDQMClient_cfg::threshold.

Referenced by doClusteringWorker().

309  {
310 
311  std::map< int, std::vector<unsigned> > hpds;
312  std::map< int, std::vector<unsigned> > rbxs;
313 
314  for(EH ih = eRecHits_.begin(); ih != eRecHits_.end(); ih++ ) {
315 
316  unsigned rhi = ih->second;
317 
318  if(! masked(rhi) ) continue;
319  // rechit was asked to be processed
320  const reco::PFRecHit& rhit = rechit( rhi, rechits);
321  //double energy = rhit.energy();
322  int layer = rhit.layer();
323  if ( layer != PFLayer::HCAL_BARREL1 &&
324  layer != PFLayer::HCAL_ENDCAP ) break;
325  // layer != PFLayer::HCAL_BARREL2) break; //BARREL2 for HO : need specific cleaning
326  HcalDetId theHcalDetId = HcalDetId(rhit.detId());
327  int ieta = theHcalDetId.ieta();
328  int iphi = theHcalDetId.iphi();
329  int ihpd = ieta < 0 ?
330  ( layer == PFLayer::HCAL_ENDCAP ? -(iphi+1)/2-100 : -iphi ) :
331  ( layer == PFLayer::HCAL_ENDCAP ? (iphi+1)/2+100 : iphi ) ;
332  hpds[ihpd].push_back(rhi);
333  int irbx = ieta < 0 ?
334  ( layer == PFLayer::HCAL_ENDCAP ? -(iphi+5)/4 - 20 : -(iphi+5)/4 ) :
335  ( layer == PFLayer::HCAL_ENDCAP ? (iphi+5)/4 + 20 : (iphi+5)/4 ) ;
336  if ( irbx == 19 ) irbx = 1;
337  else if ( irbx == -19 ) irbx = -1;
338  else if ( irbx == 39 ) irbx = 21;
339  else if ( irbx == -39 ) irbx = -21;
340  rbxs[irbx].push_back(rhi);
341  }
342 
343  // Loop on readout boxes
344  for ( std::map<int, std::vector<unsigned> >::iterator itrbx = rbxs.begin();
345  itrbx != rbxs.end(); ++itrbx ) {
346  if ( ( abs(itrbx->first)<20 && itrbx->second.size() > 30 ) ||
347  ( abs(itrbx->first)>20 && itrbx->second.size() > 30 ) ) {
348  const std::vector<unsigned>& rhits = itrbx->second;
349  double totalEta = 0.;
350  double totalEtaW = 0.;
351  double totalPhi = 0.;
352  double totalPhiW = 0.;
353  double totalEta2 = 1E-9;
354  double totalEta2W = 1E-9;
355  double totalPhi2 = 1E-9;
356  double totalPhi2W = 1E-9;
357  double totalEnergy = 0.;
358  double totalEnergy2 = 1E-9;
359  unsigned nSeeds = rhits.size();
360  unsigned nSeeds0 = rhits.size();
361  std::map< int,std::vector<unsigned> > theHPDs;
362  std::multimap< double,unsigned > theEnergies;
363  for ( unsigned jh=0; jh < rhits.size(); ++jh ) {
364  const reco::PFRecHit& hit = rechit(rhits[jh], rechits);
365  // Check if the hit is a seed
366  unsigned nN = 0;
367  bool isASeed = true;
368  const vector<unsigned>& neighbours4 = *(& hit.neighbours4());
369  for(unsigned in=0; in<neighbours4.size(); in++) {
370  const reco::PFRecHit& neighbour = rechit( neighbours4[in], rechits );
371  // one neighbour has a higher energy -> the tested rechit is not a seed
372  if( neighbour.energy() > hit.energy() ) {
373  --nSeeds;
374  --nSeeds0;
375  isASeed = false;
376  break;
377  } else {
378  if ( neighbour.energy() > 0.4 ) ++nN;
379  }
380  }
381  if ( isASeed && !nN ) --nSeeds0;
382 
383  HcalDetId theHcalDetId = HcalDetId(hit.detId());
384  // int ieta = theHcalDetId.ieta();
385  int iphi = theHcalDetId.iphi();
386  // std::cout << "Hit : " << hit.energy() << " " << ieta << " " << iphi << std::endl;
387  if ( hit.layer() == PFLayer::HCAL_BARREL1 )
388  theHPDs[iphi].push_back(rhits[jh]);
389  else
390  theHPDs[(iphi-1)/2].push_back(rhits[jh]);
391  theEnergies.insert(std::pair<double,unsigned>(hit.energy(),rhits[jh]));
392  totalEnergy += hit.energy();
393  totalPhi += fabs(hit.position().phi());
394  totalPhiW += hit.energy()*fabs(hit.position().phi());
395  totalEta += hit.position().eta();
396  totalEtaW += hit.energy()*hit.position().eta();
397  totalEnergy2 += hit.energy()*hit.energy();
398  totalPhi2 += hit.position().phi()*hit.position().phi();
399  totalPhi2W += hit.energy()*hit.position().phi()*hit.position().phi();
400  totalEta2 += hit.position().eta()*hit.position().eta();
401  totalEta2W += hit.energy()*hit.position().eta()*hit.position().eta();
402  }
403  // totalPhi /= totalEnergy;
404  totalPhi /= rhits.size();
405  totalEta /= rhits.size();
406  totalPhiW /= totalEnergy;
407  totalEtaW /= totalEnergy;
408  totalPhi2 /= rhits.size();
409  totalEta2 /= rhits.size();
410  totalPhi2W /= totalEnergy;
411  totalEta2W /= totalEnergy;
412  totalPhi2 = std::sqrt(totalPhi2 - totalPhi*totalPhi);
413  totalEta2 = std::sqrt(totalEta2 - totalEta*totalEta);
414  totalPhi2W = std::sqrt(totalPhi2W - totalPhiW*totalPhiW);
415  totalEta2W = std::sqrt(totalEta2W - totalEtaW*totalEtaW);
416  totalEnergy /= rhits.size();
417  totalEnergy2 /= rhits.size();
418  totalEnergy2 = std::sqrt(totalEnergy2 - totalEnergy*totalEnergy);
419  //if ( totalPhi2W/totalEta2W < 0.18 ) {
420  if ( nSeeds0 > 6 ) {
421  unsigned nHPD15 = 0;
422  for ( std::map<int, std::vector<unsigned> >::iterator itHPD = theHPDs.begin();
423  itHPD != theHPDs.end(); ++itHPD ) {
424  int hpdN = itHPD->first;
425  const std::vector<unsigned>& hpdHits = itHPD->second;
426  if ( ( abs(hpdN) < 100 && hpdHits.size() > 14 ) ||
427  ( abs(hpdN) > 100 && hpdHits.size() > 14 ) ) ++nHPD15;
428  }
429  if ( nHPD15 > 1 ) {
430  /*
431  std::cout << "Read out box numero " << itrbx->first
432  << " has " << itrbx->second.size() << " hits in it !"
433  << std::endl << "sigma Eta/Phi = " << totalEta2 << " " << totalPhi2 << " " << totalPhi2/totalEta2
434  << std::endl << "sigma EtaW/PhiW = " << totalEta2W << " " << totalPhi2W << " " << totalPhi2W/totalEta2W
435  << std::endl << "E = " << totalEnergy << " +/- " << totalEnergy2
436  << std::endl << "nSeeds = " << nSeeds << " " << nSeeds0
437  << std::endl;
438  for ( std::map<int, std::vector<unsigned> >::iterator itHPD = theHPDs.begin();
439  itHPD != theHPDs.end(); ++itHPD ) {
440  unsigned hpdN = itHPD->first;
441  const std::vector<unsigned>& hpdHits = itHPD->second;
442  std::cout << "HPD number " << hpdN << " contains " << hpdHits.size() << " hits" << std::endl;
443  }
444  */
445  std::multimap<double, unsigned >::iterator ntEn = theEnergies.end();
446  --ntEn;
447  unsigned nn = 0;
448  double threshold = 1.;
449  for ( std::multimap<double, unsigned >::iterator itEn = theEnergies.begin();
450  itEn != theEnergies.end(); ++itEn ) {
451  ++nn;
452  if ( nn < 5 ) {
453  mask_[itEn->second] = false;
454  } else if ( nn == 5 ) {
455  threshold = itEn->first * 5.;
456  mask_[itEn->second] = false;
457  } else {
458  if ( itEn->first < threshold ) mask_[itEn->second] = false;
459  }
460  if ( !masked(itEn->second) ) {
461  reco::PFRecHit theCleanedHit(rechit(itEn->second, rechits));
462  //theCleanedHit.setRescale(0.);
463  pfRecHitsCleaned_->push_back(theCleanedHit);
464  }
465  /*
466  if ( !masked(itEn->second) )
467  std::cout << "Hit Energies = " << itEn->first
468  << " " << ntEn->first/itEn->first << " masked " << std::endl;
469  else
470  std::cout << "Hit Energies = " << itEn->first
471  << " " << ntEn->first/itEn->first << " kept for clustering " << std::endl;
472  */
473  }
474  }
475  }
476  }
477  }
478 
479  // Loop on hpd's
480  std::map<int, std::vector<unsigned> >::iterator neighbour1;
481  std::map<int, std::vector<unsigned> >::iterator neighbour2;
482  std::map<int, std::vector<unsigned> >::iterator neighbour0;
483  std::map<int, std::vector<unsigned> >::iterator neighbour3;
484  unsigned size1 = 0;
485  unsigned size2 = 0;
486  for ( std::map<int, std::vector<unsigned> >::iterator ithpd = hpds.begin();
487  ithpd != hpds.end(); ++ithpd ) {
488 
489  const std::vector<unsigned>& rhits = ithpd->second;
490  std::multimap< double,unsigned > theEnergies;
491  double totalEnergy = 0.;
492  double totalEnergy2 = 1E-9;
493  for ( unsigned jh=0; jh < rhits.size(); ++jh ) {
494  const reco::PFRecHit& hit = rechit(rhits[jh], rechits);
495  totalEnergy += hit.energy();
496  totalEnergy2 += hit.energy()*hit.energy();
497  theEnergies.insert(std::pair<double,unsigned>(hit.energy(),rhits[jh]));
498  }
499  totalEnergy /= rhits.size();
500  totalEnergy2 /= rhits.size();
501  totalEnergy2 = std::sqrt(totalEnergy2 - totalEnergy*totalEnergy);
502 
503  if ( ithpd->first == 1 ) neighbour1 = hpds.find(72);
504  else if ( ithpd->first == -1 ) neighbour1 = hpds.find(-72);
505  else if ( ithpd->first == 101 ) neighbour1 = hpds.find(136);
506  else if ( ithpd->first == -101 ) neighbour1 = hpds.find(-136);
507  else neighbour1 = ithpd->first > 0 ? hpds.find(ithpd->first-1) : hpds.find(ithpd->first+1) ;
508 
509  if ( ithpd->first == 72 ) neighbour2 = hpds.find(1);
510  else if ( ithpd->first == -72 ) neighbour2 = hpds.find(-1);
511  else if ( ithpd->first == 136 ) neighbour2 = hpds.find(101);
512  else if ( ithpd->first == -136 ) neighbour2 = hpds.find(-101);
513  else neighbour2 = ithpd->first > 0 ? hpds.find(ithpd->first+1) : hpds.find(ithpd->first-1) ;
514 
515  if ( neighbour1 != hpds.end() ) {
516  if ( neighbour1->first == 1 ) neighbour0 = hpds.find(72);
517  else if ( neighbour1->first == -1 ) neighbour0 = hpds.find(-72);
518  else if ( neighbour1->first == 101 ) neighbour0 = hpds.find(136);
519  else if ( neighbour1->first == -101 ) neighbour0 = hpds.find(-136);
520  else neighbour0 = neighbour1->first > 0 ? hpds.find(neighbour1->first-1) : hpds.find(neighbour1->first+1) ;
521  }
522 
523  if ( neighbour2 != hpds.end() ) {
524  if ( neighbour2->first == 72 ) neighbour3 = hpds.find(1);
525  else if ( neighbour2->first == -72 ) neighbour3 = hpds.find(-1);
526  else if ( neighbour2->first == 136 ) neighbour3 = hpds.find(101);
527  else if ( neighbour2->first == -136 ) neighbour3 = hpds.find(-101);
528  else neighbour3 = neighbour2->first > 0 ? hpds.find(neighbour2->first+1) : hpds.find(neighbour2->first-1) ;
529  }
530 
531  size1 = neighbour1 != hpds.end() ? neighbour1->second.size() : 0;
532  size2 = neighbour2 != hpds.end() ? neighbour2->second.size() : 0;
533 
534  // Also treat the case of two neighbouring HPD's not in the same RBX
535  if ( size1 > 10 ) {
536  if ( ( abs(neighbour1->first) > 100 && neighbour1->second.size() > 15 ) ||
537  ( abs(neighbour1->first) < 100 && neighbour1->second.size() > 12 ) )
538  size1 = neighbour0 != hpds.end() ? neighbour0->second.size() : 0;
539  }
540  if ( size2 > 10 ) {
541  if ( ( abs(neighbour2->first) > 100 && neighbour2->second.size() > 15 ) ||
542  ( abs(neighbour2->first) < 100 && neighbour2->second.size() > 12 ) )
543  size2 = neighbour3 != hpds.end() ? neighbour3->second.size() : 0;
544  }
545 
546  if ( ( abs(ithpd->first) > 100 && ithpd->second.size() > 15 ) ||
547  ( abs(ithpd->first) < 100 && ithpd->second.size() > 12 ) )
548  if ( (float)(size1 + size2)/(float)ithpd->second.size() < 1.0 ) {
549  /*
550  std::cout << "HPD numero " << ithpd->first
551  << " has " << ithpd->second.size() << " hits in it !" << std::endl
552  << "Neighbours : " << size1 << " " << size2
553  << std::endl;
554  */
555  std::multimap<double, unsigned >::iterator ntEn = theEnergies.end();
556  --ntEn;
557  unsigned nn = 0;
558  double threshold = 1.;
559  for ( std::multimap<double, unsigned >::iterator itEn = theEnergies.begin();
560  itEn != theEnergies.end(); ++itEn ) {
561  ++nn;
562  if ( nn < 5 ) {
563  mask_[itEn->second] = false;
564  } else if ( nn == 5 ) {
565  threshold = itEn->first * 2.5;
566  mask_[itEn->second] = false;
567  } else {
568  if ( itEn->first < threshold ) mask_[itEn->second] = false;
569  }
570  if ( !masked(itEn->second) ) {
571  reco::PFRecHit theCleanedHit(rechit(itEn->second, rechits));
572  //theCleanedHit.setRescale(0.);
573  pfRecHitsCleaned_->push_back(theCleanedHit);
574  }
575  /*
576  if ( !masked(itEn->second) )
577  std::cout << "Hit Energies = " << itEn->first
578  << " " << ntEn->first/itEn->first << " masked " << std::endl;
579  else
580  std::cout << "Hit Energies = " << itEn->first
581  << " " << ntEn->first/itEn->first << " kept for clustering " << std::endl;
582  */
583  }
584  }
585  }
586 
587 }
std::multimap< double, unsigned >::iterator EH
unsigned detId() const
rechit detId
Definition: PFRecHit.h:105
const reco::PFRecHit & rechit(unsigned i, const reco::PFRecHitCollection &rechits)
std::multimap< double, unsigned, std::greater< double > > eRecHits_
indices to rechits, sorted by decreasing E (not E_T)
std::auto_ptr< std::vector< reco::PFRecHit > > pfRecHitsCleaned_
particle flow rechits cleaned
PFLayer::Layer layer() const
rechit layer
Definition: PFRecHit.h:108
bool masked(unsigned rhi) const
Particle flow rechit (rechit + geometry and topology information). See clustering algorithm in PFClus...
Definition: PFRecHit.h:31
T sqrt(T t)
Definition: SSEVec.h:48
std::vector< bool > mask_
int ieta() const
get the cell ieta
Definition: HcalDetId.h:36
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
layer definition for PFRecHit and PFCluster
Definition: PFLayer.h:21
int iphi() const
get the cell iphi
Definition: HcalDetId.h:38
const math::XYZPoint & position() const
is seed ?
Definition: PFRecHit.h:141
double energy() const
rechit energy
Definition: PFRecHit.h:111
const std::vector< unsigned > & neighbours4() const
Definition: PFRecHit.h:157
std::auto_ptr< std::vector< reco::PFCluster > >& PFClusterAlgo::clusters ( )
inline
Returns
particle flow clusters

Definition at line 220 of file PFClusterAlgo.h.

References pfClusters_.

Referenced by PFRootEventManager::clustering().

221  {return pfClusters_;}
std::auto_ptr< std::vector< reco::PFCluster > > pfClusters_
all clusters
unsigned PFClusterAlgo::color ( unsigned  rhi) const
Returns
color of the rechit

Definition at line 1766 of file PFClusterAlgo.cc.

References color_.

Referenced by cuy.plotElement::__init__(), cuy.superimposeElement::__init__(), cuy.graphElement::__init__(), DisplayManager::loadGRecHits(), and paint().

1766  {
1767 
1768  if(rhi>=color_.size() ) { // rhi >= 0, since rhi is unsigned
1769  string err = "PFClusterAlgo::color : out of range";
1770  throw std::out_of_range(err);
1771  }
1772 
1773  return color_[rhi];
1774 }
std::vector< unsigned > color_
color, for all rechits
reco::PFRecHitRef PFClusterAlgo::createRecHitRef ( const reco::PFRecHitCollection rechits,
unsigned  rhi 
)
private

create a reference to a rechit. in case rechitsHandle_.isValid(), this reference is permanent.

Definition at line 1801 of file PFClusterAlgo.cc.

References edm::HandleBase::isValid(), and rechitsHandle_.

Referenced by buildPFClusters().

1802  {
1803 
1804  if( rechitsHandle_.isValid() ) {
1805  return reco::PFRecHitRef( rechitsHandle_, rhi );
1806  }
1807  else {
1808  return reco::PFRecHitRef( &rechits, rhi );
1809  }
1810 }
PFRecHitHandle rechitsHandle_
edm::Ref< PFRecHitCollection > PFRecHitRef
persistent reference to PFRecHit objects
Definition: PFRecHitFwd.h:15
bool isValid() const
Definition: HandleBase.h:76
std::pair< double, double > PFClusterAlgo::dCrack ( double  phi,
double  eta 
)
private

distance to a crack in the ECAL barrel in eta and phi direction

Definition at line 1851 of file PFClusterAlgo.cc.

References gather_cfg::cout, alignCSCRings::e, i, M_PI, bookConverter::min, and pi.

Referenced by findSeeds().

1851  {
1852 
1853  static const double pi= M_PI;// 3.14159265358979323846;
1854 
1855  //Location of the 18 phi-cracks
1856  static std::vector<double> cPhi;
1857  if(cPhi.size()==0)
1858  {
1859  cPhi.resize(18,0);
1860  cPhi[0]=2.97025;
1861  for(unsigned i=1;i<=17;++i) cPhi[i]=cPhi[0]-2*i*pi/18;
1862  }
1863 
1864  //Shift of this location if eta<0
1865  static double delta_cPhi=0.00638;
1866 
1867  double defi; //the result
1868 
1869  //the location is shifted
1870  if(eta<0) phi +=delta_cPhi;
1871 
1872  if (phi>=-pi && phi<=pi){
1873 
1874  //the problem of the extrema
1875  if (phi<cPhi[17] || phi>=cPhi[0]){
1876  if (phi<0) phi+= 2*pi;
1877  defi = std::min(fabs(phi -cPhi[0]),fabs(phi-cPhi[17]-2*pi));
1878  }
1879 
1880  //between these extrema...
1881  else{
1882  bool OK = false;
1883  unsigned i=16;
1884  while(!OK){
1885  if (phi<cPhi[i]){
1886  defi=std::min(fabs(phi-cPhi[i+1]),fabs(phi-cPhi[i]));
1887  OK=true;
1888  }
1889  else i-=1;
1890  }
1891  }
1892  }
1893  else{
1894  defi=0.; //if there is a problem, we assum that we are in a crack
1895  std::cout<<"Problem in dminphi"<<std::endl;
1896  }
1897  //if(eta<0) defi=-defi; //because of the disymetry
1898 
1899  static std::vector<double> cEta;
1900  if ( cEta.size() == 0 ) {
1901  cEta.push_back(0.0);
1902  cEta.push_back(4.44747e-01) ; cEta.push_back(-4.44747e-01) ;
1903  cEta.push_back(7.92824e-01) ; cEta.push_back(-7.92824e-01) ;
1904  cEta.push_back(1.14090e+00) ; cEta.push_back(-1.14090e+00) ;
1905  cEta.push_back(1.47464e+00) ; cEta.push_back(-1.47464e+00) ;
1906  }
1907  double deta = 999.; // the other result
1908 
1909  for ( unsigned ieta=0; ieta<cEta.size(); ++ieta ) {
1910  deta = std::min(deta,fabs(eta-cEta[ieta]));
1911  }
1912 
1913  defi /= 0.0175;
1914  deta /= 0.0175;
1915  return std::pair<double,double>(defi,deta);
1916 
1917 }
int i
Definition: DBlmapReader.cc:9
pair< int, edm::FunctionWithDict > OK
Definition: findMethod.cc:70
T eta() const
#define M_PI
Definition: BFit3D.cc:3
tuple cout
Definition: gather_cfg.py:121
double pi
Definition: DDAxes.h:10
void PFClusterAlgo::doClustering ( const reco::PFRecHitCollection rechits)

perform clustering

Definition at line 114 of file PFClusterAlgo.cc.

References edm::HandleBase::clear(), doClusteringWorker(), kNotDefined, mask_, rechitsHandle_, sortedRecHits_, funct::true, and which_pos_calc_.

Referenced by PFRootEventManager::clustering().

114  {
115 
116  // using rechits without a Handle, clear to avoid a stale member
118  if( which_pos_calc_ != kNotDefined ) {
120  sortedRecHits_->sort();
121  }
122 
123  // clear rechits mask
124  mask_.clear();
125  mask_.resize( rechits.size(), true );
126 
127  // perform clustering
129 
130 }
PFRecHitHandle rechitsHandle_
void doClusteringWorker(const reco::PFRecHitCollection &rechits)
perform clustering
std::vector< bool > mask_
PositionCalcType which_pos_calc_
parameter for position calculation
std::auto_ptr< SortedPFRecHitCollection > sortedRecHits_
edm::SortedCollection< reco::PFRecHit > SortedPFRecHitCollection
Definition: PFClusterAlgo.h:58
void PFClusterAlgo::doClustering ( const reco::PFRecHitCollection rechits,
const std::vector< bool > &  mask 
)

Definition at line 132 of file PFClusterAlgo.cc.

References edm::HandleBase::clear(), doClusteringWorker(), kNotDefined, mask_, rechitsHandle_, sortedRecHits_, funct::true, and which_pos_calc_.

132  {
133  // using rechits without a Handle, clear to avoid a stale member
134 
136  if( which_pos_calc_ != kNotDefined ) {
138  sortedRecHits_->sort();
139  }
140 
141  // use the specified mask, unless it doesn't match with the rechits
142  mask_.clear();
143 
144  if (mask.size() == rechits.size()) {
145  mask_.insert( mask_.end(), mask.begin(), mask.end() );
146  } else {
147  edm::LogError("PClusterAlgo::doClustering") << "map size should be " << rechits.size() << ". Will be reinitialized.";
148  mask_.resize( rechits.size(), true );
149  }
150 
151  // perform clustering
153 
154 }
PFRecHitHandle rechitsHandle_
void doClusteringWorker(const reco::PFRecHitCollection &rechits)
perform clustering
std::vector< bool > mask_
PositionCalcType which_pos_calc_
parameter for position calculation
std::auto_ptr< SortedPFRecHitCollection > sortedRecHits_
edm::SortedCollection< reco::PFRecHit > SortedPFRecHitCollection
Definition: PFClusterAlgo.h:58
void PFClusterAlgo::doClustering ( const PFRecHitHandle rechitsHandle)

perform clustering in full framework

Definition at line 68 of file PFClusterAlgo.cc.

References doClusteringWorker(), kNotDefined, mask_, HI_PhotonSkim_cff::rechits, rechitsHandle_, sortedRecHits_, funct::true, and which_pos_calc_.

68  {
69  const reco::PFRecHitCollection& rechits = * rechitsHandle;
70 
71  // cache the Handle to the rechits
72  rechitsHandle_ = rechitsHandle;
73  if( which_pos_calc_ != kNotDefined ) {
75  sortedRecHits_->sort();
76  }
77 
78  // clear rechits mask
79  mask_.clear();
80  mask_.resize( rechits.size(), true );
81 
82  // perform clustering
83  doClusteringWorker( rechits );
84 }
PFRecHitHandle rechitsHandle_
std::vector< PFRecHit > PFRecHitCollection
collection of PFRecHit objects
Definition: PFRecHitFwd.h:9
void doClusteringWorker(const reco::PFRecHitCollection &rechits)
perform clustering
std::vector< bool > mask_
PositionCalcType which_pos_calc_
parameter for position calculation
std::auto_ptr< SortedPFRecHitCollection > sortedRecHits_
edm::SortedCollection< reco::PFRecHit > SortedPFRecHitCollection
Definition: PFClusterAlgo.h:58
void PFClusterAlgo::doClustering ( const PFRecHitHandle rechitsHandle,
const std::vector< bool > &  mask 
)

Definition at line 86 of file PFClusterAlgo.cc.

References doClusteringWorker(), kNotDefined, mask_, HI_PhotonSkim_cff::rechits, rechitsHandle_, sortedRecHits_, funct::true, and which_pos_calc_.

87  {
88 
89  const reco::PFRecHitCollection& rechits = * rechitsHandle;
90 
91  // cache the Handle to the rechits
92  rechitsHandle_ = rechitsHandle;
93  if( which_pos_calc_ != kNotDefined ) {
95  sortedRecHits_->sort();
96  }
97 
98 
99  // use the specified mask, unless it doesn't match with the rechits
100  mask_.clear();
101  if (mask.size() == rechits.size()) {
102  mask_.insert( mask_.end(), mask.begin(), mask.end() );
103  } else {
104  edm::LogError("PClusterAlgo::doClustering") << "map size should be " << rechits.size() << ". Will be reinitialized.";
105  mask_.resize( rechits.size(), true );
106  }
107 
108  // perform clustering
109 
110  doClusteringWorker( rechits );
111 
112 }
PFRecHitHandle rechitsHandle_
std::vector< PFRecHit > PFRecHitCollection
collection of PFRecHit objects
Definition: PFRecHitFwd.h:9
void doClusteringWorker(const reco::PFRecHitCollection &rechits)
perform clustering
std::vector< bool > mask_
PositionCalcType which_pos_calc_
parameter for position calculation
std::auto_ptr< SortedPFRecHitCollection > sortedRecHits_
edm::SortedCollection< reco::PFRecHit > SortedPFRecHitCollection
Definition: PFClusterAlgo.h:58
void PFClusterAlgo::doClusteringWorker ( const reco::PFRecHitCollection rechits)
private

perform clustering

Definition at line 156 of file PFClusterAlgo.cc.

References buildPFClusters(), buildTopoClusters(), cleanRBXAndHPD(), cleanRBXandHPDs_, color_, relval_parameters_module::energy, eRecHits_, edm::false, findSeeds(), i, pfClusters_, pfRecHitsCleaned_, rechit(), seedStates_, topoClusters_, UNKNOWN, and usedInTopo_.

Referenced by doClustering().

156  {
157 
158 
159  if ( pfClusters_.get() )
160  pfClusters_->clear();
161  else
162  pfClusters_.reset( new std::vector<reco::PFCluster> );
163 
164 
165  if ( pfRecHitsCleaned_.get() )
166  pfRecHitsCleaned_->clear();
167  else
168  pfRecHitsCleaned_.reset( new std::vector<reco::PFRecHit> );
169 
170 
171  eRecHits_.clear();
172 
173  for ( unsigned i = 0; i < rechits.size(); i++ ){
174 
175  eRecHits_.insert( make_pair( rechit(i, rechits).energy(), i) );
176  }
177 
178  color_.clear();
179  color_.resize( rechits.size(), 0 );
180 
181  seedStates_.clear();
182  seedStates_.resize( rechits.size(), UNKNOWN );
183 
184  usedInTopo_.clear();
185  usedInTopo_.resize( rechits.size(), false );
186 
188 
189  // look for seeds.
190 
191  findSeeds( rechits );
192 
193  // build topological clusters around seeds
195 
196  // look for PFClusters inside each topological cluster (one per seed)
197 
198 
199  // int ix=0;
200  // for (reco::PFRecHitCollection::const_iterator cand =rechits.begin(); cand<rechits.end(); cand++){
201  // cout <<ix++<<" "<< cand->layer()<<endl;
202  // }
203 
204 
205  for(unsigned i=0; i<topoClusters_.size(); i++) {
206 
207  const std::vector< unsigned >& topocluster = topoClusters_[i];
208 
209  buildPFClusters( topocluster, rechits );
210 
211  }
212 
213 }
void cleanRBXAndHPD(const reco::PFRecHitCollection &rechits)
Clean HCAL readout box noise and HPD discharge.
int i
Definition: DBlmapReader.cc:9
std::vector< unsigned > color_
color, for all rechits
const reco::PFRecHit & rechit(unsigned i, const reco::PFRecHitCollection &rechits)
std::multimap< double, unsigned, std::greater< double > > eRecHits_
indices to rechits, sorted by decreasing E (not E_T)
std::vector< SeedState > seedStates_
seed state, for all rechits
std::vector< std::vector< unsigned > > topoClusters_
sets of cells having one common side, and energy over threshold
std::auto_ptr< std::vector< reco::PFRecHit > > pfRecHitsCleaned_
particle flow rechits cleaned
std::vector< bool > usedInTopo_
used in topo cluster? for all rechits
bool cleanRBXandHPDs_
option to clean HCAL RBX&#39;s and HPD&#39;s
void buildTopoClusters(const reco::PFRecHitCollection &rechits)
build topoclusters around seeds
void findSeeds(const reco::PFRecHitCollection &rechits)
look for seeds
std::auto_ptr< std::vector< reco::PFCluster > > pfClusters_
all clusters
volatile std::atomic< bool > shutdown_flag false
void buildPFClusters(const std::vector< unsigned > &cluster, const reco::PFRecHitCollection &rechits)
build PFClusters from a topocluster
void PFClusterAlgo::enableDebugging ( bool  debug)
inline

set hits on which clustering will be performed

enable/disable debugging

Definition at line 76 of file PFClusterAlgo.h.

References debug, and debug_.

Referenced by PFRootEventManager::readOptions().

76 { debug_ = debug;}
bool debug_
debugging on/off
#define debug
Definition: HDRShower.cc:19
void PFClusterAlgo::findSeeds ( const reco::PFRecHitCollection rechits)
private

look for seeds

Definition at line 590 of file PFClusterAlgo.cc.

References funct::abs(), dtNoiseDBValidation_cfg::cerr, CLEAN, CLEAN_S4S1, CLEAN_THRESH, gather_cfg::cout, dCrack(), debug_, DOUBLESPIKE_S6S2, DOUBLESPIKE_THRESH, PFLayer::ECAL_BARREL, PFLayer::ECAL_ENDCAP, reco::PFRecHit::energy(), reco::PFRecHit::energyUp(), eRecHits_, eta(), file_, hBNeighbour, PFLayer::HCAL_BARREL1, PFLayer::HCAL_BARREL2, PFLayer::HCAL_ENDCAP, hENeighbour, PFLayer::HF_EM, PFLayer::HF_HAD, recoMuon::in, reco::PFRecHit::layer(), mask_, masked(), bookConverter::min, reco::PFRecHit::neighbours4(), reco::PFRecHit::neighbours8(), nNeighbours_, NO, paint(), parameter(), pfRecHitsCleaned_, phi, reco::PFRecHit::position(), reco::PFRecHit::positionREP(), PFLayer::PS1, PFLayer::PS2, reco::PFRecHit::pt2(), rechit(), SEED, SEED_PT_THRESH, SEED_THRESH, seeds_, seedStates_, and YES.

Referenced by doClusteringWorker().

590  {
591 
592  seeds_.clear();
593 
594  // should replace this by the message logger.
595 #ifdef PFLOW_DEBUG
596  if(debug_)
597  cout<<"PFClusterAlgo::findSeeds : start"<<endl;
598 #endif
599 
600 
601  // An empty list of neighbours
602  const vector<unsigned> noNeighbours(0, static_cast<unsigned>(0));
603 
604  // loop on rechits (sorted by decreasing energy - not E_T)
605 
606  for(EH ih = eRecHits_.begin(); ih != eRecHits_.end(); ih++ ) {
607 
608  unsigned rhi = ih->second;
609 
610  if(! masked(rhi) ) continue;
611  // rechit was asked to be processed
612 
613  double rhenergy = ih->first;
614  const reco::PFRecHit& wannaBeSeed = rechit(rhi, rechits);
615 
616  if( seedStates_[rhi] == NO ) continue;
617  // this hit was already tested, and is not a seed
618 
619  // determine seed energy threshold depending on the detector
620  int layer = wannaBeSeed.layer();
621  //for HO Ring0 and 1/2 boundary
622 
623  int iring = 0;
624  if (layer==PFLayer::HCAL_BARREL2 && abs(wannaBeSeed.positionREP().Eta())>0.34) iring= 1;
625 
626  double seedThresh = parameter( SEED_THRESH,
627  static_cast<PFLayer::Layer>(layer), 0, iring );
628 
629  double seedPtThresh = parameter( SEED_PT_THRESH,
630  static_cast<PFLayer::Layer>(layer), 0., iring );
631 
632  double cleanThresh = parameter( CLEAN_THRESH,
633  static_cast<PFLayer::Layer>(layer), 0, iring );
634 
635  double minS4S1_a = parameter( CLEAN_S4S1,
636  static_cast<PFLayer::Layer>(layer), 0, iring );
637 
638  double minS4S1_b = parameter( CLEAN_S4S1,
639  static_cast<PFLayer::Layer>(layer), 1, iring );
640 
641  double doubleSpikeThresh = parameter( DOUBLESPIKE_THRESH,
642  static_cast<PFLayer::Layer>(layer), 0, iring );
643 
644  double doubleSpikeS6S2 = parameter( DOUBLESPIKE_S6S2,
645  static_cast<PFLayer::Layer>(layer), 0, iring );
646 
647 #ifdef PFLOW_DEBUG
648  if(debug_)
649  cout<<"layer:"<<layer<<" seedThresh:"<<seedThresh<<endl;
650 #endif
651 
652 
653  if( rhenergy < seedThresh || (seedPtThresh>0. && wannaBeSeed.pt2() < seedPtThresh*seedPtThresh )) {
654  seedStates_[rhi] = NO;
655  continue;
656  }
657 
658  // Find the cell unused neighbours
659  const vector<unsigned>* nbp;
660  double tighterE = 1.0;
661  double tighterF = 1.0;
662 
663  switch ( layer ) {
664  case PFLayer::ECAL_BARREL:
665  case PFLayer::ECAL_ENDCAP:
669  tighterE = 2.0;
670  tighterF = 3.0;
671  case PFLayer::HF_EM:
672  case PFLayer::HF_HAD:
673  if( nNeighbours_ == 4 ) {
674  nbp = & wannaBeSeed.neighbours4();
675  }
676  else if( nNeighbours_ == 8 ) {
677  nbp = & wannaBeSeed.neighbours8();
678  }
679  else if( nNeighbours_ == 0 ) {
680  nbp = & noNeighbours;
681  // Allows for no clustering at all: all rechits are clusters.
682  // Useful for HF
683  }
684  else {
685  cerr<<"you're not allowed to set n neighbours to "
686  <<nNeighbours_<<endl;
687  assert(0);
688  }
689  break;
690  case PFLayer::PS1:
691  case PFLayer::PS2:
692  nbp = & wannaBeSeed.neighbours4();
693  break;
694 
695  default:
696  cerr<<"CellsEF::PhotonSeeds : unknown layer "<<layer<<endl;
697  assert(0);
698  }
699 
700  const vector<unsigned>& neighbours = *nbp;
701 
702 
703  // Select as a seed if all neighbours have a smaller energy
704 
705  seedStates_[rhi] = YES;
706  for(unsigned in=0; in<neighbours.size(); in++) {
707 
708  unsigned rhj = neighbours[in];
709  // Ignore neighbours already masked
710  if ( !masked(rhj) ) continue;
711  const reco::PFRecHit& neighbour = rechit( rhj, rechits );
712 
713  // one neighbour has a higher energy -> the tested rechit is not a seed
714  if( neighbour.energy() > wannaBeSeed.energy() ) {
715  seedStates_[rhi] = NO;
716  break;
717  }
718  }
719 
720  // Cleaning : check energetic, isolated seeds, likely to come from erratic noise.
721  if ( file_ || wannaBeSeed.energy() > cleanThresh ) {
722 
723  const vector<unsigned>& neighbours4 = *(& wannaBeSeed.neighbours4());
724  // Determine the fraction of surrounding energy
725  double surroundingEnergy = wannaBeSeed.energyUp();
726  double neighbourEnergy = 0.;
727  double layerEnergy = 0.;
728  for(unsigned in4=0; in4<neighbours4.size(); in4++) {
729  unsigned rhj = neighbours4[in4];
730  // Ignore neighbours already masked
731  if ( !masked(rhj) ) continue;
732  const reco::PFRecHit& neighbour = rechit( rhj, rechits );
733  surroundingEnergy += neighbour.energy() + neighbour.energyUp();
734  neighbourEnergy += neighbour.energy() + neighbour.energyUp();
735  layerEnergy += neighbour.energy();
736  }
737  // Fraction 0 is the balance between EM and HAD layer for this tower
738  // double fraction0 = layer == PFLayer::HF_EM || layer == PFLayer::HF_HAD ?
739  // wannaBeSeed.energyUp()/wannaBeSeed.energy() : 1.;
740  // Fraction 1 is the balance between the hit and its neighbours from both layers
741  double fraction1 = surroundingEnergy/wannaBeSeed.energy();
742  // Fraction 2 is the balance between the tower and the tower neighbours
743  // double fraction2 = neighbourEnergy/(wannaBeSeed.energy()+wannaBeSeed.energyUp());
744  // Fraction 3 is the balance between the hits and the hits neighbours in the same layer.
745  // double fraction3 = layerEnergy/(wannaBeSeed.energy());
746  // Mask the seed and the hit if energetic/isolated rechit
747  // if ( fraction0 < minS4S1 || fraction1 < minS4S1 || fraction2 < minS4S1 || fraction3 < minS4S1 ) {
748  // if ( fraction1 < minS4S1 || ( wannaBeSeed.energy() > 1.5*cleanThresh && fraction0 + fraction3 < minS4S1 ) ) {
749 
750  if ( file_ ) {
751  if ( layer == PFLayer::ECAL_BARREL || layer == PFLayer::HCAL_BARREL1 || layer == PFLayer::HCAL_BARREL2) { //BARREL2 for HO
752  /*
753  double eta = wannaBeSeed.position().eta();
754  double phi = wannaBeSeed.position().phi();
755  std::pair<double,double> dcr = dCrack(phi,eta);
756  double dcrmin = std::min(dcr.first, dcr.second);
757  if ( dcrmin > 1. )
758  */
759  hBNeighbour->Fill(1./wannaBeSeed.energy(), fraction1);
760  } else if ( fabs(wannaBeSeed.position().eta()) < 5.0 ) {
761  if ( layer == PFLayer::ECAL_ENDCAP || layer == PFLayer::HCAL_ENDCAP ) {
762  if ( wannaBeSeed.energy() > 1000 ) {
763  if ( fabs(wannaBeSeed.position().eta()) < 2.8 )
764  hENeighbour->Fill(1./wannaBeSeed.energy(), fraction1);
765  }
766  } else {
767  hENeighbour->Fill(1./wannaBeSeed.energy(), fraction1);
768  }
769  }
770  }
771 
772  if ( wannaBeSeed.energy() > cleanThresh ) {
773  double f1Cut = minS4S1_a * log10(wannaBeSeed.energy()) + minS4S1_b;
774  if ( fraction1 < f1Cut ) {
775  // Double the energy cleaning threshold when close to the ECAL/HCAL - HF transition
776  double eta = wannaBeSeed.position().eta();
777  double phi = wannaBeSeed.position().phi();
778  std::pair<double,double> dcr = dCrack(phi,eta);
779  double dcrmin = layer == PFLayer::ECAL_BARREL ? std::min(dcr.first, dcr.second) : dcr.second;
780  eta = fabs(eta);
781  if ( eta < 5.0 && // No cleaning for the HF border
782  ( ( eta < 2.85 && dcrmin > 1. ) ||
783  ( rhenergy > tighterE*cleanThresh &&
784  fraction1 < f1Cut/tighterF ) ) // Tighter cleaning for various cracks
785  ) {
786  seedStates_[rhi] = CLEAN;
787  mask_[rhi] = false;
788  reco::PFRecHit theCleanedHit(wannaBeSeed);
789  //theCleanedHit.setRescale(0.);
790  pfRecHitsCleaned_->push_back(theCleanedHit);
791  /*
792  std::cout << "A seed with E/pT/eta/phi = " << wannaBeSeed.energy() << " " << wannaBeSeed.energyUp()
793  << " " << sqrt(wannaBeSeed.pt2()) << " " << wannaBeSeed.position().eta() << " " << phi
794  << " and with surrounding fractions = " << fraction0 << " " << fraction1
795  << " " << fraction2 << " " << fraction3
796  << " in layer " << layer
797  << " had been cleaned " << std::endl
798  << " Distances to cracks : " << dcr.first << " " << dcr.second << std::endl
799  << "(Cuts were : " << cleanThresh << " and " << f1Cut << ")" << std::endl;
800  */
801  }
802  }
803  }
804  }
805 
806  // Clean double spikes
807  if ( mask_[rhi] && wannaBeSeed.energy() > doubleSpikeThresh ) {
808  // Determine energy surrounding the seed and the most energetic neighbour
809  double surroundingEnergyi = 0.;
810  double enmax = -999.;
811  unsigned mostEnergeticNeighbour = 0;
812  const vector<unsigned>& neighbours4i = *(& wannaBeSeed.neighbours4());
813  for(unsigned in4=0; in4<neighbours4i.size(); in4++) {
814  unsigned rhj = neighbours4i[in4];
815  if ( !masked(rhj) ) continue;
816  const reco::PFRecHit& neighbouri = rechit( rhj, rechits );
817  surroundingEnergyi += neighbouri.energy();
818  if ( neighbouri.energy() > enmax ) {
819  enmax = neighbouri.energy();
820  mostEnergeticNeighbour = rhj;
821  }
822  }
823  // Is there an energetic neighbour ?
824  if ( enmax > 0. ) {
825  unsigned rhj = mostEnergeticNeighbour;
826  const reco::PFRecHit& neighbouri = rechit( rhj, rechits );
827  double surroundingEnergyj = 0.;
828  //if ( mask_[rhj] && neighbouri.energy() > doubleSpikeThresh ) {
829  // Determine energy surrounding the energetic neighbour
830  const vector<unsigned>& neighbours4j = *(& neighbouri.neighbours4());
831  for(unsigned jn4=0; jn4<neighbours4j.size(); jn4++) {
832  unsigned rhk = neighbours4j[jn4];
833  const reco::PFRecHit& neighbourj = rechit( rhk, rechits );
834  surroundingEnergyj += neighbourj.energy();
835  }
836  // The energy surrounding the double spike candidate
837  double surroundingEnergyFraction =
838  (surroundingEnergyi+surroundingEnergyj) / (wannaBeSeed.energy()+neighbouri.energy()) - 1.;
839  if ( surroundingEnergyFraction < doubleSpikeS6S2 ) {
840  double eta = wannaBeSeed.position().eta();
841  double phi = wannaBeSeed.position().phi();
842  std::pair<double,double> dcr = dCrack(phi,eta);
843  double dcrmin = layer == PFLayer::ECAL_BARREL ? std::min(dcr.first, dcr.second) : dcr.second;
844  eta = fabs(eta);
845  if ( ( eta < 5.0 && dcrmin > 1. ) ||
846  ( wannaBeSeed.energy() > tighterE*doubleSpikeThresh &&
847  surroundingEnergyFraction < doubleSpikeS6S2/tighterF ) ) {
848  /*
849  std::cout << "Double spike cleaned : Energies = " << wannaBeSeed.energy()
850  << " " << neighbouri.energy()
851  << " surrounded by only " << surroundingEnergyFraction*100.
852  << "% of the two spike energies at eta/phi/distance to closest crack = "
853  << eta << " " << phi << " " << dcrmin
854  << std::endl;
855  */
856  // mask the seed
857  seedStates_[rhi] = CLEAN;
858  mask_[rhi] = false;
859  reco::PFRecHit theCleanedSeed(wannaBeSeed);
860  pfRecHitsCleaned_->push_back(theCleanedSeed);
861  // mask the neighbour
862  seedStates_[rhj] = CLEAN;
863  mask_[rhj] = false;
864  reco::PFRecHit theCleanedNeighbour(wannaBeSeed);
865  pfRecHitsCleaned_->push_back(neighbouri);
866  }
867  }
868  } else {
869  /*
870  std::cout << "PFClusterAlgo : Double spike search : An isolated cell should have been killed " << std::endl
871  << "but is going through the double spike search!" << std::endl;
872  */
873  }
874  }
875 
876  if ( seedStates_[rhi] == YES ) {
877 
878  // seeds_ contains the indices of all seeds.
879  seeds_.push_back( rhi );
880 
881  // marking the rechit
882  paint(rhi, SEED);
883 
884  // then all neighbours cannot be seeds and are flagged as such
885  for(unsigned in=0; in<neighbours.size(); in++) {
886  seedStates_[ neighbours[in] ] = NO;
887  }
888  }
889 
890  }
891 
892 #ifdef PFLOW_DEBUG
893  if(debug_)
894  cout<<"PFClusterAlgo::findSeeds : done"<<endl;
895 #endif
896 }
std::multimap< double, unsigned >::iterator EH
double pt2() const
rechit momentum transverse to the beam, squared.
Definition: PFRecHit.h:123
std::pair< double, double > dCrack(double phi, double eta)
distance to a crack in the ECAL barrel in eta and phi direction
const reco::PFRecHit & rechit(unsigned i, const reco::PFRecHitCollection &rechits)
TH2F * hENeighbour
std::multimap< double, unsigned, std::greater< double > > eRecHits_
indices to rechits, sorted by decreasing E (not E_T)
T eta() const
std::vector< SeedState > seedStates_
seed state, for all rechits
std::auto_ptr< std::vector< reco::PFRecHit > > pfRecHitsCleaned_
particle flow rechits cleaned
double energyUp() const
For HF hits: rechit energy (and neighbour&#39;s) in the other HF layer.
Definition: PFRecHit.h:120
PFLayer::Layer layer() const
rechit layer
Definition: PFRecHit.h:108
U second(std::pair< T, U > const &p)
bool masked(unsigned rhi) const
std::vector< unsigned > seeds_
vector of indices for seeds.
TH2F * hBNeighbour
Particle flow rechit (rechit + geometry and topology information). See clustering algorithm in PFClus...
Definition: PFRecHit.h:31
const std::vector< unsigned > & neighbours8() const
Definition: PFRecHit.h:160
std::vector< bool > mask_
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
const math::XYZPoint & position() const
is seed ?
Definition: PFRecHit.h:141
bool debug_
debugging on/off
int nNeighbours_
number of neighbours
void paint(unsigned rhi, unsigned color=1)
paint a rechit with a color.
double parameter(Parameter paramtype, PFLayer::Layer layer, unsigned iCoeff=0, int iring0=0) const
double energy() const
rechit energy
Definition: PFRecHit.h:111
tuple cout
Definition: gather_cfg.py:121
const std::vector< unsigned > & neighbours4() const
Definition: PFRecHit.h:157
const REPPoint & positionREP() const
rechit cell centre rho, eta, phi. call calculatePositionREP before !
Definition: PFRecHit.cc:121
Definition: DDAxes.h:10
bool PFClusterAlgo::isSeed ( unsigned  rhi) const
Returns
seed flag (not seed state) for rechit with index rhi

Definition at line 1778 of file PFClusterAlgo.cc.

References seedStates_.

Referenced by buildPFClusters(), calculateClusterPosition(), and PFRootEventManager::printRecHits().

1778  {
1779 
1780  if(rhi>=seedStates_.size() ) { // rhi >= 0, since rhi is unsigned
1781  string err = "PFClusterAlgo::isSeed : out of range";
1782  throw std::out_of_range(err);
1783  }
1784 
1785  return seedStates_[rhi]>0 ? true : false;
1786 }
std::vector< SeedState > seedStates_
seed state, for all rechits
bool PFClusterAlgo::masked ( unsigned  rhi) const
Returns
mask flag for the rechit

Definition at line 1755 of file PFClusterAlgo.cc.

References mask_.

Referenced by buildTopoCluster(), buildTopoClusters(), cleanRBXAndHPD(), and findSeeds().

1755  {
1756 
1757  if(rhi>=mask_.size() ) { // rhi >= 0, since rhi is unsigned
1758  string err = "PFClusterAlgo::masked : out of range";
1759  throw std::out_of_range(err);
1760  }
1761 
1762  return mask_[rhi];
1763 }
std::vector< bool > mask_
int PFClusterAlgo::nNeighbours ( ) const
inline

get number of neighbours for

Definition at line 183 of file PFClusterAlgo.h.

References nNeighbours_.

183 { return nNeighbours_;}
int nNeighbours_
number of neighbours
void PFClusterAlgo::paint ( unsigned  rhi,
unsigned  color = 1 
)
private

paint a rechit with a color.

Definition at line 1789 of file PFClusterAlgo.cc.

References color(), and color_.

Referenced by buildPFClusters(), and findSeeds().

1789  {
1790 
1791  if(rhi>=color_.size() ) { // rhi >= 0, since rhi is unsigned
1792  string err = "PFClusterAlgo::color : out of range";
1793  throw std::out_of_range(err);
1794  }
1795 
1796  color_[rhi] = color;
1797 }
std::vector< unsigned > color_
color, for all rechits
unsigned color(unsigned rhi) const
double PFClusterAlgo::parameter ( Parameter  paramtype,
PFLayer::Layer  layer,
unsigned  iCoeff = 0,
int  iring0 = 0 
) const
Returns
the value of a parameter of a given type, see enum Parameter, in a given layer.

Definition at line 216 of file PFClusterAlgo.cc.

References dtNoiseDBValidation_cfg::cerr, CLEAN_S4S1, CLEAN_THRESH, DOUBLESPIKE_S6S2, DOUBLESPIKE_THRESH, PFLayer::ECAL_BARREL, PFLayer::ECAL_ENDCAP, PFLayer::HCAL_BARREL1, PFLayer::HCAL_BARREL2, PFLayer::HCAL_ENDCAP, PFLayer::HF_EM, PFLayer::HF_HAD, minS4S1Barrel_, minS4S1Endcap_, minS6S2DoubleSpikeBarrel_, minS6S2DoubleSpikeEndcap_, PFLayer::PS1, PFLayer::PS2, PT_THRESH, SEED_PT_THRESH, SEED_THRESH, THRESH, threshBarrel_, threshCleanBarrel_, threshCleanEndcap_, threshDoubleSpikeBarrel_, threshDoubleSpikeEndcap_, threshEndcap_, threshPtBarrel_, threshPtEndcap_, threshPtSeedBarrel_, threshPtSeedEndcap_, threshSeedBarrel_, threshSeedEndcap_, and relativeConstraints::value.

Referenced by buildTopoCluster(), and findSeeds().

218  {
219 
220 
221  double value = 0;
222 
223  switch( layer ) {
224 
227  case PFLayer::HCAL_BARREL2: // HO.
228 
229  switch(paramtype) {
230  case THRESH:
231  value = (iring==0) ? threshBarrel_ : threshEndcap_; //For HO Ring0 and others
232  break;
233  case SEED_THRESH:
234  value = (iring==0) ? threshSeedBarrel_ : threshSeedEndcap_;
235  break;
236  case PT_THRESH:
237  value = (iring==0) ? threshPtBarrel_ : threshPtEndcap_;
238  break;
239  case SEED_PT_THRESH:
240  value = (iring==0) ? threshPtSeedBarrel_ : threshPtSeedEndcap_;
241  break;
242  case CLEAN_THRESH:
243  value = threshCleanBarrel_;
244  break;
245  case CLEAN_S4S1:
246  value = minS4S1Barrel_[iCoeff];
247  break;
248  case DOUBLESPIKE_THRESH:
249  value = threshDoubleSpikeBarrel_;
250  break;
251  case DOUBLESPIKE_S6S2:
253  break;
254  default:
255  cerr<<"PFClusterAlgo::parameter : unknown parameter type "
256  <<paramtype<<endl;
257  assert(0);
258  }
259  break;
262  case PFLayer::PS1:
263  case PFLayer::PS2:
264  case PFLayer::HF_EM:
265  case PFLayer::HF_HAD:
266  // and no particle flow in VFCAL
267  switch(paramtype) {
268  case THRESH:
269  value = threshEndcap_;
270  break;
271  case SEED_THRESH:
272  value = threshSeedEndcap_;
273  break;
274  case PT_THRESH:
275  value = threshPtEndcap_;
276  break;
277  case SEED_PT_THRESH:
278  value = threshPtSeedEndcap_;
279  break;
280  case CLEAN_THRESH:
281  value = threshCleanEndcap_;
282  break;
283  case CLEAN_S4S1:
284  value = minS4S1Endcap_[iCoeff];
285  break;
286  case DOUBLESPIKE_THRESH:
287  value = threshDoubleSpikeEndcap_;
288  break;
289  case DOUBLESPIKE_S6S2:
291  break;
292  default:
293  cerr<<"PFClusterAlgo::parameter : unknown parameter type "
294  <<paramtype<<endl;
295  assert(0);
296  }
297  break;
298  default:
299  cerr<<"PFClusterAlgo::parameter : unknown layer "<<layer<<endl;
300  assert(0);
301  break;
302  }
303 
304  return value;
305 }
double threshEndcap_
endcap threshold
double threshDoubleSpikeBarrel_
Barrel double-spike cleaning.
double threshDoubleSpikeEndcap_
Endcap double-spike cleaning.
double threshPtEndcap_
std::vector< double > minS4S1Barrel_
double threshPtSeedEndcap_
double threshPtSeedBarrel_
double threshSeedBarrel_
barrel seed threshold
std::vector< double > minS4S1Endcap_
double threshCleanBarrel_
Barrel cleaning threshold and S4/S1 smallest fractiom.
double minS6S2DoubleSpikeEndcap_
double threshBarrel_
barrel threshold
double threshSeedEndcap_
endcap seed threshold
double threshPtBarrel_
double minS6S2DoubleSpikeBarrel_
double threshCleanEndcap_
Endcap cleaning threshold and S4/S1 smallest fractiom.
int PFClusterAlgo::posCalcNCrystal ( ) const
inline

get number of crystals for position calculation (-1 all,5, or 9)

Definition at line 189 of file PFClusterAlgo.h.

References posCalcNCrystal_.

Referenced by buildPFClusters().

189 {return posCalcNCrystal_;}
int posCalcNCrystal_
number of crystals for position calculation
double PFClusterAlgo::posCalcP1 ( ) const
inline

get p1 for position calculation

Definition at line 186 of file PFClusterAlgo.h.

References posCalcP1_.

186 { return posCalcP1_; }
const reco::PFRecHit & PFClusterAlgo::rechit ( unsigned  i,
const reco::PFRecHitCollection rechits 
)


Returns
hit with index i. performs a test on the vector range

Definition at line 1742 of file PFClusterAlgo.cc.

References i.

Referenced by buildPFClusters(), buildTopoCluster(), cleanRBXAndHPD(), doClusteringWorker(), and findSeeds().

1743  {
1744 
1745  if(i >= rechits.size() ) { // i >= 0, since i is unsigned
1746  string err = "PFClusterAlgo::rechit : out of range";
1747  throw std::out_of_range(err);
1748  }
1749 
1750  return rechits[i];
1751 }
int i
Definition: DBlmapReader.cc:9
std::auto_ptr< std::vector< reco::PFRecHit > >& PFClusterAlgo::rechitsCleaned ( )
inline
Returns
cleaned rechits

Definition at line 224 of file PFClusterAlgo.h.

References pfRecHitsCleaned_.

225  {return pfRecHitsCleaned_;}
std::auto_ptr< std::vector< reco::PFRecHit > > pfRecHitsCleaned_
particle flow rechits cleaned
void PFClusterAlgo::setCleanRBXandHPDs ( bool  cleanRBXandHPDs)
inline

Activate cleaning of HCAL RBX's and HPD's.

Definition at line 164 of file PFClusterAlgo.h.

References cleanRBXandHPDs_.

Referenced by PFRootEventManager::readOptions().

164 { cleanRBXandHPDs_ = cleanRBXandHPDs; }
bool cleanRBXandHPDs_
option to clean HCAL RBX&#39;s and HPD&#39;s
void PFClusterAlgo::setEBGeom ( const CaloSubdetectorGeometry esh)
inline

Definition at line 138 of file PFClusterAlgo.h.

References eb_geom.

138  {
139  eb_geom = esh;
140  }
const CaloSubdetectorGeometry * eb_geom
void PFClusterAlgo::setEEGeom ( const CaloSubdetectorGeometry esh)
inline

Definition at line 141 of file PFClusterAlgo.h.

References ee_geom.

141  {
142  ee_geom = esh;
143  }
const CaloSubdetectorGeometry * ee_geom
void PFClusterAlgo::setEGammaPosCalc ( const edm::ParameterSet conf)
inline

Definition at line 135 of file PFClusterAlgo.h.

References eg_pos_calc.

135  {
136  eg_pos_calc.reset(new PositionCalc(conf));
137  }
std::auto_ptr< PositionCalc > eg_pos_calc
void PFClusterAlgo::setHistos ( TFile *  file,
TH2F *  hB,
TH2F *  hE 
)
inline

set endcap clean threshold

Definition at line 126 of file PFClusterAlgo.h.

References mergeVDriftHistosByStation::file, file_, hBNeighbour, and hENeighbour.

Referenced by PFRootEventManager::readOptions().

void PFClusterAlgo::setNNeighbours ( int  n)
inline

set number of neighbours for

Definition at line 129 of file PFClusterAlgo.h.

References n, and nNeighbours_.

Referenced by PFRootEventManager::readOptions().

129 { nNeighbours_ = n;}
int nNeighbours_
number of neighbours
void PFClusterAlgo::setPosCalcNCrystal ( int  n)
inline

set number of crystals for position calculation (-1 all,5, or 9)

Definition at line 155 of file PFClusterAlgo.h.

References n, and posCalcNCrystal_.

Referenced by PFRootEventManager::readOptions().

155 { posCalcNCrystal_ = n;}
int posCalcNCrystal_
number of crystals for position calculation
void PFClusterAlgo::setPosCalcP1 ( double  p1)
inline

set p1 for position calculation

Definition at line 152 of file PFClusterAlgo.h.

References p1, and posCalcP1_.

Referenced by PFRootEventManager::readOptions().

152 { posCalcP1_ = p1; }
double p1[4]
Definition: TauolaWrapper.h:89
void PFClusterAlgo::setPosCalcW0 ( double  w0)
inline

Definition at line 149 of file PFClusterAlgo.h.

References param_W0_.

149 { param_W0_ = w0; }
void PFClusterAlgo::setPositionCalcType ( const PositionCalcType t)
inline

Definition at line 132 of file PFClusterAlgo.h.

References lumiQTWidget::t, and which_pos_calc_.

132 {which_pos_calc_ = t;}
PositionCalcType which_pos_calc_
parameter for position calculation
void PFClusterAlgo::setPreshowerGeom ( const CaloSubdetectorGeometry esh)
inline

Definition at line 144 of file PFClusterAlgo.h.

References preshower_geom.

144  {
145  preshower_geom = esh;
146  }
const CaloSubdetectorGeometry * preshower_geom
void PFClusterAlgo::setS4S1CleanBarrel ( const std::vector< double > &  coeffs)
inline

Definition at line 103 of file PFClusterAlgo.h.

References minS4S1Barrel_.

Referenced by PFRootEventManager::readOptions().

103 {minS4S1Barrel_ = coeffs;}
std::vector< double > minS4S1Barrel_
void PFClusterAlgo::setS4S1CleanEndcap ( const std::vector< double > &  coeffs)
inline

Definition at line 119 of file PFClusterAlgo.h.

References minS4S1Endcap_.

Referenced by PFRootEventManager::readOptions().

119 {minS4S1Endcap_ = coeffs;}
std::vector< double > minS4S1Endcap_
void PFClusterAlgo::setS6S2DoubleSpikeBarrel ( double  cut)
inline
void PFClusterAlgo::setS6S2DoubleSpikeEndcap ( double  cut)
inline
void PFClusterAlgo::setShowerSigma ( double  sigma)
inline

set shower sigma for

Definition at line 158 of file PFClusterAlgo.h.

References showerSigma_.

Referenced by PFRootEventManager::readOptions().

158 { showerSigma_ = sigma;}
double showerSigma_
sigma of shower (cm)
void PFClusterAlgo::setThreshBarrel ( double  thresh)
inline

setters -------------------------------------------------——

set barrel threshold

Definition at line 94 of file PFClusterAlgo.h.

References GOODCOLL_filter_cfg::thresh, and threshBarrel_.

Referenced by PFRootEventManager::readOptions().

double threshBarrel_
barrel threshold
void PFClusterAlgo::setThreshCleanBarrel ( double  thresh)
inline

set barrel clean threshold

Definition at line 102 of file PFClusterAlgo.h.

References GOODCOLL_filter_cfg::thresh, and threshCleanBarrel_.

Referenced by PFRootEventManager::readOptions().

double threshCleanBarrel_
Barrel cleaning threshold and S4/S1 smallest fractiom.
void PFClusterAlgo::setThreshCleanEndcap ( double  thresh)
inline

set endcap clean threshold

Definition at line 118 of file PFClusterAlgo.h.

References GOODCOLL_filter_cfg::thresh, and threshCleanEndcap_.

Referenced by PFRootEventManager::readOptions().

double threshCleanEndcap_
Endcap cleaning threshold and S4/S1 smallest fractiom.
void PFClusterAlgo::setThreshDoubleSpikeBarrel ( double  thresh)
inline

set endcap thresholds for double spike cleaning

Definition at line 106 of file PFClusterAlgo.h.

References GOODCOLL_filter_cfg::thresh, and threshDoubleSpikeBarrel_.

Referenced by PFRootEventManager::readOptions().

double threshDoubleSpikeBarrel_
Barrel double-spike cleaning.
void PFClusterAlgo::setThreshDoubleSpikeEndcap ( double  thresh)
inline

set endcap thresholds for double spike cleaning

Definition at line 122 of file PFClusterAlgo.h.

References GOODCOLL_filter_cfg::thresh, and threshDoubleSpikeEndcap_.

Referenced by PFRootEventManager::readOptions().

double threshDoubleSpikeEndcap_
Endcap double-spike cleaning.
void PFClusterAlgo::setThreshEndcap ( double  thresh)
inline

set endcap threshold

Definition at line 110 of file PFClusterAlgo.h.

References GOODCOLL_filter_cfg::thresh, and threshEndcap_.

Referenced by PFRootEventManager::readOptions().

double threshEndcap_
endcap threshold
void PFClusterAlgo::setThreshPtBarrel ( double  thresh)
inline
void PFClusterAlgo::setThreshPtEndcap ( double  thresh)
inline
void PFClusterAlgo::setThreshPtSeedBarrel ( double  thresh)
inline
void PFClusterAlgo::setThreshPtSeedEndcap ( double  thresh)
inline
void PFClusterAlgo::setThreshSeedBarrel ( double  thresh)
inline

set barrel seed threshold

Definition at line 98 of file PFClusterAlgo.h.

References GOODCOLL_filter_cfg::thresh, and threshSeedBarrel_.

Referenced by PFRootEventManager::readOptions().

double threshSeedBarrel_
barrel seed threshold
void PFClusterAlgo::setThreshSeedEndcap ( double  thresh)
inline

set endcap seed threshold

Definition at line 114 of file PFClusterAlgo.h.

References GOODCOLL_filter_cfg::thresh, and threshSeedEndcap_.

Referenced by PFRootEventManager::readOptions().

double threshSeedEndcap_
endcap seed threshold
void PFClusterAlgo::setUseCornerCells ( bool  usecornercells)
inline

activate use of cells with a common corner to build topo-clusters

Definition at line 161 of file PFClusterAlgo.h.

References useCornerCells_.

Referenced by PFRootEventManager::readOptions().

161 { useCornerCells_ = usecornercells;}
bool useCornerCells_
option to use cells with a common corner to build topo-clusters
double PFClusterAlgo::showerSigma ( ) const
inline

get shower sigma

Definition at line 192 of file PFClusterAlgo.h.

References showerSigma_.

192 { return showerSigma_ ;}
double showerSigma_
sigma of shower (cm)
double PFClusterAlgo::threshBarrel ( ) const
inline

getters -------------------------------------------------——

get barrel threshold

Definition at line 169 of file PFClusterAlgo.h.

References threshBarrel_.

Referenced by DisplayManager::createGRecHit().

169 {return threshBarrel_;}
double threshBarrel_
barrel threshold
double PFClusterAlgo::threshEndcap ( ) const
inline

get endcap threshold

Definition at line 176 of file PFClusterAlgo.h.

References threshEndcap_.

Referenced by DisplayManager::createGRecHit().

176 {return threshEndcap_;}
double threshEndcap_
endcap threshold
double PFClusterAlgo::threshSeedBarrel ( ) const
inline

get barrel seed threshold

Definition at line 172 of file PFClusterAlgo.h.

References threshSeedBarrel_.

172 {return threshSeedBarrel_;}
double threshSeedBarrel_
barrel seed threshold
double PFClusterAlgo::threshSeedEndcap ( ) const
inline

get endcap seed threshold

Definition at line 179 of file PFClusterAlgo.h.

References threshSeedEndcap_.

179 {return threshSeedEndcap_;}
double threshSeedEndcap_
endcap seed threshold
void PFClusterAlgo::write ( void  )

write histos

Definition at line 57 of file PFClusterAlgo.cc.

References gather_cfg::cout, and file_.

Referenced by PFRootEventManager::write().

57  {
58 
59  if ( file_ ) {
60  file_->Write();
61  cout << "Benchmark output written to file " << file_->GetName() << endl;
62  file_->Close();
63  }
64 
65 }
tuple cout
Definition: gather_cfg.py:121

Friends And Related Function Documentation

std::ostream& operator<< ( std::ostream &  out,
const PFClusterAlgo algo 
)
friend

Member Data Documentation

bool PFClusterAlgo::cleanRBXandHPDs_
private

option to clean HCAL RBX's and HPD's

Definition at line 387 of file PFClusterAlgo.h.

Referenced by doClusteringWorker(), and setCleanRBXandHPDs().

std::vector< unsigned > PFClusterAlgo::color_
private

color, for all rechits

Definition at line 312 of file PFClusterAlgo.h.

Referenced by color(), doClusteringWorker(), and paint().

bool PFClusterAlgo::debug_
private

debugging on/off

Definition at line 390 of file PFClusterAlgo.h.

Referenced by buildPFClusters(), buildTopoCluster(), buildTopoClusters(), enableDebugging(), and findSeeds().

const CaloSubdetectorGeometry* PFClusterAlgo::eb_geom
private

Definition at line 378 of file PFClusterAlgo.h.

Referenced by calculateClusterPosition(), and setEBGeom().

const CaloSubdetectorGeometry * PFClusterAlgo::ee_geom
private

Definition at line 378 of file PFClusterAlgo.h.

Referenced by calculateClusterPosition(), and setEEGeom().

std::auto_ptr<PositionCalc> PFClusterAlgo::eg_pos_calc
private

Definition at line 377 of file PFClusterAlgo.h.

Referenced by calculateClusterPosition(), and setEGammaPosCalc().

std::multimap<double, unsigned, std::greater<double> > PFClusterAlgo::eRecHits_
private

indices to rechits, sorted by decreasing E (not E_T)

Definition at line 305 of file PFClusterAlgo.h.

Referenced by cleanRBXAndHPD(), doClusteringWorker(), and findSeeds().

TFile* PFClusterAlgo::file_
private

Definition at line 399 of file PFClusterAlgo.h.

Referenced by findSeeds(), PFClusterAlgo(), setHistos(), and write().

TH2F* PFClusterAlgo::hBNeighbour
private

Definition at line 397 of file PFClusterAlgo.h.

Referenced by findSeeds(), PFClusterAlgo(), and setHistos().

TH2F* PFClusterAlgo::hENeighbour
private

Definition at line 398 of file PFClusterAlgo.h.

Referenced by findSeeds(), PFClusterAlgo(), and setHistos().

std::set<unsigned> PFClusterAlgo::idUsedRecHits_
private

ids of rechits used in seed search

Definition at line 302 of file PFClusterAlgo.h.

std::vector< bool > PFClusterAlgo::mask_
private

mask, for all rechits. only masked rechits will be clustered by default, all rechits are masked. see setMask function.

Definition at line 309 of file PFClusterAlgo.h.

Referenced by cleanRBXAndHPD(), doClustering(), findSeeds(), and masked().

std::vector<double> PFClusterAlgo::minS4S1Barrel_
private

Definition at line 353 of file PFClusterAlgo.h.

Referenced by operator<<(), parameter(), and setS4S1CleanBarrel().

std::vector<double> PFClusterAlgo::minS4S1Endcap_
private

Definition at line 361 of file PFClusterAlgo.h.

Referenced by operator<<(), parameter(), and setS4S1CleanEndcap().

double PFClusterAlgo::minS6S2DoubleSpikeBarrel_
private

Definition at line 357 of file PFClusterAlgo.h.

Referenced by parameter(), and setS6S2DoubleSpikeBarrel().

double PFClusterAlgo::minS6S2DoubleSpikeEndcap_
private

Definition at line 365 of file PFClusterAlgo.h.

Referenced by parameter(), and setS6S2DoubleSpikeEndcap().

int PFClusterAlgo::nNeighbours_
private

number of neighbours

Definition at line 368 of file PFClusterAlgo.h.

Referenced by findSeeds(), nNeighbours(), operator<<(), and setNNeighbours().

double PFClusterAlgo::param_W0_
private

Definition at line 376 of file PFClusterAlgo.h.

Referenced by calculateClusterPosition(), and setPosCalcW0().

std::auto_ptr< std::vector<reco::PFCluster> > PFClusterAlgo::pfClusters_
private

all clusters

particle flow clusters

Definition at line 330 of file PFClusterAlgo.h.

Referenced by buildPFClusters(), clusters(), doClusteringWorker(), and operator<<().

std::auto_ptr< std::vector<reco::PFRecHit> > PFClusterAlgo::pfRecHitsCleaned_
private

particle flow rechits cleaned

Definition at line 333 of file PFClusterAlgo.h.

Referenced by cleanRBXAndHPD(), doClusteringWorker(), findSeeds(), and rechitsCleaned().

int PFClusterAlgo::posCalcNCrystal_
private

number of crystals for position calculation

Definition at line 371 of file PFClusterAlgo.h.

Referenced by buildPFClusters(), calculateClusterPosition(), operator<<(), posCalcNCrystal(), and setPosCalcNCrystal().

double PFClusterAlgo::posCalcP1_
private

Definition at line 375 of file PFClusterAlgo.h.

Referenced by calculateClusterPosition(), operator<<(), posCalcP1(), and setPosCalcP1().

const CaloSubdetectorGeometry * PFClusterAlgo::preshower_geom
private

Definition at line 378 of file PFClusterAlgo.h.

Referenced by calculateClusterPosition(), and setPreshowerGeom().

unsigned PFClusterAlgo::prodNum_ = 1
staticprivate

product number

Definition at line 394 of file PFClusterAlgo.h.

PFRecHitHandle PFClusterAlgo::rechitsHandle_
private

Definition at line 297 of file PFClusterAlgo.h.

Referenced by createRecHitRef(), and doClustering().

std::vector< unsigned > PFClusterAlgo::seeds_
private

vector of indices for seeds.

Definition at line 321 of file PFClusterAlgo.h.

Referenced by buildTopoClusters(), and findSeeds().

std::vector< SeedState > PFClusterAlgo::seedStates_
private

seed state, for all rechits

Definition at line 315 of file PFClusterAlgo.h.

Referenced by buildPFClusters(), doClusteringWorker(), findSeeds(), and isSeed().

double PFClusterAlgo::showerSigma_
private

sigma of shower (cm)

Definition at line 381 of file PFClusterAlgo.h.

Referenced by buildPFClusters(), operator<<(), setShowerSigma(), and showerSigma().

std::auto_ptr<SortedPFRecHitCollection> PFClusterAlgo::sortedRecHits_
private

Definition at line 299 of file PFClusterAlgo.h.

Referenced by calculateClusterPosition(), and doClustering().

double PFClusterAlgo::threshBarrel_
private

barrel threshold

Definition at line 336 of file PFClusterAlgo.h.

Referenced by calculateClusterPosition(), operator<<(), parameter(), setThreshBarrel(), and threshBarrel().

double PFClusterAlgo::threshCleanBarrel_
private

Barrel cleaning threshold and S4/S1 smallest fractiom.

Definition at line 352 of file PFClusterAlgo.h.

Referenced by operator<<(), parameter(), and setThreshCleanBarrel().

double PFClusterAlgo::threshCleanEndcap_
private

Endcap cleaning threshold and S4/S1 smallest fractiom.

Definition at line 360 of file PFClusterAlgo.h.

Referenced by operator<<(), parameter(), and setThreshCleanEndcap().

double PFClusterAlgo::threshDoubleSpikeBarrel_
private

Barrel double-spike cleaning.

Definition at line 356 of file PFClusterAlgo.h.

Referenced by parameter(), and setThreshDoubleSpikeBarrel().

double PFClusterAlgo::threshDoubleSpikeEndcap_
private

Endcap double-spike cleaning.

Definition at line 364 of file PFClusterAlgo.h.

Referenced by parameter(), and setThreshDoubleSpikeEndcap().

double PFClusterAlgo::threshEndcap_
private

endcap threshold

Definition at line 344 of file PFClusterAlgo.h.

Referenced by calculateClusterPosition(), operator<<(), parameter(), setThreshEndcap(), and threshEndcap().

double PFClusterAlgo::threshPtBarrel_
private

Definition at line 337 of file PFClusterAlgo.h.

Referenced by operator<<(), parameter(), and setThreshPtBarrel().

double PFClusterAlgo::threshPtEndcap_
private

Definition at line 345 of file PFClusterAlgo.h.

Referenced by operator<<(), parameter(), and setThreshPtEndcap().

double PFClusterAlgo::threshPtSeedBarrel_
private

Definition at line 341 of file PFClusterAlgo.h.

Referenced by operator<<(), parameter(), and setThreshPtSeedBarrel().

double PFClusterAlgo::threshPtSeedEndcap_
private

Definition at line 349 of file PFClusterAlgo.h.

Referenced by operator<<(), parameter(), and setThreshPtSeedEndcap().

double PFClusterAlgo::threshSeedBarrel_
private

barrel seed threshold

Definition at line 340 of file PFClusterAlgo.h.

Referenced by operator<<(), parameter(), setThreshSeedBarrel(), and threshSeedBarrel().

double PFClusterAlgo::threshSeedEndcap_
private

endcap seed threshold

Definition at line 348 of file PFClusterAlgo.h.

Referenced by operator<<(), parameter(), setThreshSeedEndcap(), and threshSeedEndcap().

std::vector< std::vector< unsigned > > PFClusterAlgo::topoClusters_
private

sets of cells having one common side, and energy over threshold

Definition at line 324 of file PFClusterAlgo.h.

Referenced by buildTopoClusters(), and doClusteringWorker().

bool PFClusterAlgo::useCornerCells_
private

option to use cells with a common corner to build topo-clusters

Definition at line 384 of file PFClusterAlgo.h.

Referenced by buildTopoCluster(), operator<<(), and setUseCornerCells().

std::vector< bool > PFClusterAlgo::usedInTopo_
private

used in topo cluster? for all rechits

Definition at line 318 of file PFClusterAlgo.h.

Referenced by buildTopoCluster(), buildTopoClusters(), and doClusteringWorker().

PositionCalcType PFClusterAlgo::which_pos_calc_
private

parameter for position calculation

Definition at line 374 of file PFClusterAlgo.h.

Referenced by calculateClusterPosition(), doClustering(), and setPositionCalcType().