Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
RawClusterBuilderTopo.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file RawClusterBuilderTopo.cc
1 #include "RawClusterBuilderTopo.h"
2 
3 #include <calobase/RawCluster.h>
4 #include <calobase/RawClusterContainer.h>
5 #include <calobase/RawClusterv1.h>
6 #include <calobase/RawTowerDefs.h> // for encode_towerid, Calorime...
7 #include <calobase/RawTowerGeom.h>
8 #include <calobase/RawTowerGeomContainer.h>
9 #include <calobase/TowerInfo.h> // for TowerInfo
10 #include <calobase/TowerInfoContainer.h> // for TowerInfoContainer
11 
12 #include <fun4all/Fun4AllReturnCodes.h>
13 #include <fun4all/SubsysReco.h>
14 
15 #include <phool/PHCompositeNode.h>
16 #include <phool/PHIODataNode.h>
17 #include <phool/PHNode.h>
18 #include <phool/PHNodeIterator.h>
19 #include <phool/PHObject.h>
20 #include <phool/getClass.h>
21 #include <phool/phool.h>
22 
23 #include <algorithm>
24 #include <cmath>
25 #include <cstdlib> // for abs
26 #include <exception>
27 #include <iostream>
28 #include <iterator> // for begin, end
29 #include <list>
30 #include <memory> // for allocator_traits<>::valu...
31 #include <stdexcept>
32 #include <utility>
33 #include <vector>
34 
35 bool sort_by_pair_second(const std::pair<int, float> &a, const std::pair<int, float> &b)
36 {
37  return (a.second > b.second);
38 }
39 
40 int RawClusterBuilderTopo::RawClusterBuilderTopo_constants_EMCal_eta_start_given_IHCal[24] =
41  {2, 6, 10, 14, 18, 22, 26, 30, 33, 37, 41, 44,
42  48, 52, 55, 59, 63, 66, 70, 74, 78, 82, 86, 90};
43 
44 int RawClusterBuilderTopo::RawClusterBuilderTopo_constants_EMCal_eta_end_given_IHCal[24] =
45  {5, 9, 13, 17, 21, 25, 29, 32, 36, 40, 43, 47,
46  51, 54, 58, 62, 65, 69, 73, 77, 81, 85, 89, 93};
47 
48 int RawClusterBuilderTopo::RawClusterBuilderTopo_constants_IHCal_eta_given_EMCal[96] = {
49  -1, -1, 0, 0, 0, 0, 1, 1, 1, 1, 2, 2,
50  2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 5, 5,
51  5, 5, 6, 6, 6, 6, 7, 7, 7, 8, 8, 8,
52  8, 9, 9, 9, 9, 10, 10, 10, 11, 11, 11, 11,
53  12, 12, 12, 12, 13, 13, 13, 14, 14, 14, 14, 15,
54  15, 15, 15, 16, 16, 16, 17, 17, 17, 17, 18, 18,
55  18, 18, 19, 19, 19, 19, 20, 20, 20, 20, 21, 21,
56  21, 21, 22, 22, 22, 22, 23, 23, 23, 23, -1, -1};
57 
58 float RawClusterBuilderTopo::calculate_dR(float eta1, float eta2, float phi1, float phi2)
59 {
60  float deta = eta1 - eta2;
61  float dphi = phi1 - phi2;
62  while (dphi > M_PI)
63  {
64  dphi -= 2 * M_PI;
65  }
66  while (dphi < -M_PI)
67  {
68  dphi += 2 * M_PI;
69  }
70  return sqrt(pow(deta, 2) + pow(dphi, 2));
71 }
72 
73 std::vector<int> RawClusterBuilderTopo::get_adjacent_towers_by_ID(int ID)
74 {
75  int this_layer = get_ilayer_from_ID(ID);
76  int this_eta = get_ieta_from_ID(ID);
77  int this_phi = get_iphi_from_ID(ID);
78 
79  std::vector<int> adjacent_towers;
80 
81  // for both IHCal and OHCal, add adjacent layers in the HCal
82  if (this_layer == 0 || this_layer == 1)
83  {
84  for (int delta_layer = 0; delta_layer <= 1; delta_layer++)
85  {
86  for (int delta_eta = -1; delta_eta <= 1; delta_eta++)
87  {
88  for (int delta_phi = -1; delta_phi <= 1; delta_phi++)
89  {
90  if (delta_layer == 0 && delta_eta == 0 && delta_phi == 0)
91  {
92  continue; // this is the same tower
93  }
94 
95  int test_eta = this_eta + delta_eta;
96  if (test_eta < 0 || test_eta >= _HCAL_NETA)
97  {
98  continue;
99  } // ignore if at the (eta) edge of calorimeter
100 
101  int test_layer = (this_layer + delta_layer) % 2; // wrap around in layer
102  int test_phi = (this_phi + delta_phi + _HCAL_NPHI) % _HCAL_NPHI; // wrap around in phi (add 64 to avoid -1)
103 
104  // disallow "corner" adjacency (diagonal in eta/phi plane and in different layer) if this option not enabled
105  if (!_allow_corner_neighbor && delta_layer == 1 && abs(delta_phi) == 1 && abs(delta_eta) == 1)
106  {
107  if (Verbosity() > 20)
108  {
109  std::cout << "RawClusterBuilderTopo::get_adjacent_towers_by_ID : corner growth not allowed " << std::endl;
110  }
111  continue;
112  }
113 
114  // add to list of adjacent towers
115  adjacent_towers.push_back(get_ID(test_layer, test_eta, test_phi));
116  }
117  }
118  }
119  }
120 
121  // for IHCal only, also add 4x4 group of EMCal towers
122  if (this_layer == 0 && _enable_EMCal)
123  {
124  int EMCal_phi_start = get_first_matching_EMCal_phi_from_IHCal(this_phi);
125  int EMCal_eta_start = RawClusterBuilderTopo_constants_EMCal_eta_start_given_IHCal[this_eta];
126  int EMCal_eta_end = RawClusterBuilderTopo_constants_EMCal_eta_end_given_IHCal[this_eta];
127 
128  for (int new_eta = EMCal_eta_start; new_eta <= EMCal_eta_end; new_eta++)
129  {
130  for (int delta_phi = 0; delta_phi < 4; delta_phi++)
131  {
132  int new_phi = (EMCal_phi_start + delta_phi + _EMCAL_NPHI) % _EMCAL_NPHI;
133 
134  int EMCal_tower = get_ID(2, new_eta, new_phi);
135  if (Verbosity() > 20)
136  {
137  std::cout << "RawClusterBuilderTopo::get_adjacent_towers_by_ID : HCal tower with eta / phi = " << this_eta << " / " << this_phi << ", adding EMCal tower with eta / phi = " << new_eta << " / " << new_phi << std::endl;
138  }
139  adjacent_towers.push_back(EMCal_tower);
140  }
141  }
142  }
143 
144  // for EMCal, add adjacent EMCal towers and (usually) one IHCal tower
145  if (this_layer == 2)
146  {
147  // EMCal towers first
148  for (int delta_eta = -1; delta_eta <= 1; delta_eta++)
149  {
150  for (int delta_phi = -1; delta_phi <= 1; delta_phi++)
151  {
152  if (delta_eta == 0 && delta_phi == 0)
153  {
154  continue; // this is the same tower
155  }
156 
157  int test_eta = this_eta + delta_eta;
158  if (test_eta < 0 || test_eta >= _EMCAL_NETA)
159  {
160  continue;
161  } // ignore if at the (eta) edge of calorimeter
162 
163  int test_phi = (this_phi + delta_phi + _EMCAL_NPHI) % _EMCAL_NPHI; // wrap around in phi (add 256 to avoid -1)
164 
165  // add to list of adjacent towers
166  adjacent_towers.push_back(get_ID(this_layer, test_eta, test_phi));
167  }
168  }
169 
170  // now add IHCal towers
171  if (_enable_HCal)
172  {
173  int HCal_eta = RawClusterBuilderTopo_constants_IHCal_eta_given_EMCal[this_eta];
174  int HCal_phi = get_matching_HCal_phi_from_EMCal(this_phi);
175 
176  if (HCal_eta >= 0)
177  {
178  int IHCal_tower = get_ID(0, HCal_eta, HCal_phi);
179  if (Verbosity() > 20)
180  {
181  std::cout << "RawClusterBuilderTopo::get_adjacent_towers_by_ID : EMCal tower with eta / phi = " << this_eta << " / " << this_phi << ", adding IHCal tower with eta / phi = " << HCal_eta << " / " << HCal_phi << std::endl;
182  }
183  adjacent_towers.push_back(IHCal_tower);
184  }
185  else
186  {
187  if (Verbosity() > 20)
188  {
189  std::cout << "RawClusterBuilderTopo::get_adjacent_towers_by_ID : EMCal tower with eta / phi = " << this_eta << " / " << this_phi << ", does not have matching IHCal due to large eta " << std::endl;
190  }
191  }
192  }
193  }
194 
195  return adjacent_towers;
196 }
197 
198 void RawClusterBuilderTopo::export_single_cluster(const std::vector<int> &original_towers)
199 {
200  if (Verbosity() > 2)
201  {
202  std::cout << "RawClusterBuilderTopo::export_single_cluster called " << std::endl;
203  }
204 
205  std::map<int, std::pair<int, int> > tower_ownership;
206  for (const int &original_tower : original_towers)
207  {
208  tower_ownership[original_tower] = std::pair<int, int>(0, -1); // all towers owned by cluster 0
209  }
210  export_clusters(original_towers, tower_ownership, 1, std::vector<float>(), std::vector<float>(), std::vector<float>());
211 
212  return;
213 }
214 
215 void RawClusterBuilderTopo::export_clusters(const std::vector<int> &original_towers, std::map<int, std::pair<int, int> > tower_ownership, unsigned int n_clusters, std::vector<float> pseudocluster_sumE, std::vector<float> pseudocluster_eta, std::vector<float> pseudocluster_phi)
216 {
217  if (n_clusters != 1) // if we didn't just pass down from export_single_cluster
218  {
219  if (Verbosity() > 2)
220  {
221  std::cout << "RawClusterBuilderTopo::export_clusters called on an initial cluster with " << n_clusters << " final clusters " << std::endl;
222  }
223  }
224  // build a RawCluster for output
225  std::vector<RawCluster *> clusters;
226  std::vector<float> clusters_E;
227  std::vector<float> clusters_x;
228  std::vector<float> clusters_y;
229  std::vector<float> clusters_z;
230 
231  for (unsigned int pc = 0; pc < n_clusters; pc++)
232  {
233  clusters.push_back(new RawClusterv1());
234  clusters_E.push_back(0);
235  clusters_x.push_back(0);
236  clusters_y.push_back(0);
237  clusters_z.push_back(0);
238  }
239 
240  for (int original_tower : original_towers)
241  {
242  int this_ID = original_tower;
243  std::pair<int, int> the_pair = tower_ownership[this_ID];
244 
245  if (Verbosity() > 5)
246  {
247  std::cout << "RawClusterBuilderTopo::export_clusters -> assigning tower " << original_tower << " with ownership ( " << the_pair.first << ", " << the_pair.second << " ) " << std::endl;
248  }
249  int this_ieta = get_ieta_from_ID(this_ID);
250  int this_iphi = get_iphi_from_ID(this_ID);
251  int this_layer = get_ilayer_from_ID(this_ID);
252  float this_E = get_E_from_ID(this_ID);
253 
254  int this_key = 0;
255  if (this_layer == 2)
256  {
257  this_key = _EMTOWERMAP_KEY_ETA_PHI[this_ieta][this_iphi];
258  }
259  else
260  {
261  this_key = _TOWERMAP_KEY_LAYER_ETA_PHI[this_layer][this_ieta][this_iphi];
262  }
263 
264  RawTowerGeom *tower_geom = _geom_containers[this_layer]->get_tower_geometry(this_key);
265 
266  if (the_pair.second == -1)
267  {
268  // assigned only to one cluster, easy
269  clusters[the_pair.first]->addTower(this_key, this_E);
270  clusters_E[the_pair.first] = clusters_E[the_pair.first] + this_E;
271  clusters_x[the_pair.first] = clusters_x[the_pair.first] + this_E * tower_geom->get_center_x();
272  clusters_y[the_pair.first] = clusters_y[the_pair.first] + this_E * tower_geom->get_center_y();
273  clusters_z[the_pair.first] = clusters_z[the_pair.first] + this_E * tower_geom->get_center_z();
274 
275  if (Verbosity() > 5)
276  {
277  std::cout << " -> tower ID " << this_ID << " fully assigned to pseudocluster " << the_pair.first << std::endl;
278  }
279  }
280  else
281  {
282  // assigned to two clusters! get energy sharing fraction ...
283  float dR1 = calculate_dR(tower_geom->get_eta(), pseudocluster_eta[the_pair.first], tower_geom->get_phi(), pseudocluster_phi[the_pair.first]) / _R_shower;
284  float dR2 = calculate_dR(tower_geom->get_eta(), pseudocluster_eta[the_pair.second], tower_geom->get_phi(), pseudocluster_phi[the_pair.second]) / _R_shower;
285  float r = std::exp(dR1 - dR2);
286  float frac1 = pseudocluster_sumE[the_pair.first] / (pseudocluster_sumE[the_pair.first] + r * pseudocluster_sumE[the_pair.second]);
287 
288  if (Verbosity() > 5)
289  {
290  std::cout << " tower ID " << this_ID << " has dR1 = " << dR1 << " to pseudocluster " << the_pair.first << " , and dR2 = " << dR2 << " to pseudocluster " << the_pair.second << ", so frac1 = " << frac1 << std::endl;
291  }
292  clusters[the_pair.first]->addTower(this_key, this_E * frac1);
293  clusters_E[the_pair.first] = clusters_E[the_pair.first] + this_E * frac1;
294  clusters_x[the_pair.first] = clusters_x[the_pair.first] + this_E * tower_geom->get_center_x() * frac1;
295  clusters_y[the_pair.first] = clusters_y[the_pair.first] + this_E * tower_geom->get_center_y() * frac1;
296  clusters_z[the_pair.first] = clusters_z[the_pair.first] + this_E * tower_geom->get_center_z() * frac1;
297 
298  clusters[the_pair.second]->addTower(this_key, this_E * (1 - frac1));
299  clusters_E[the_pair.second] = clusters_E[the_pair.second] + this_E * (1 - frac1);
300  clusters_x[the_pair.second] = clusters_x[the_pair.second] + this_E * tower_geom->get_center_x() * (1 - frac1);
301  clusters_y[the_pair.second] = clusters_y[the_pair.second] + this_E * tower_geom->get_center_y() * (1 - frac1);
302  clusters_z[the_pair.second] = clusters_z[the_pair.second] + this_E * tower_geom->get_center_z() * (1 - frac1);
303  }
304  }
305 
306  // iterate through and add to official container
307 
308  for (unsigned int cl = 0; cl < n_clusters; cl++)
309  {
310  clusters[cl]->set_energy(clusters_E[cl]);
311 
312  float mean_x = clusters_x[cl] / clusters_E[cl];
313  float mean_y = clusters_y[cl] / clusters_E[cl];
314  float mean_z = clusters_z[cl] / clusters_E[cl];
315 
316  clusters[cl]->set_r(std::sqrt(mean_y * mean_y + mean_x * mean_x));
317  clusters[cl]->set_phi(std::atan2(mean_y, mean_x));
318  clusters[cl]->set_z(mean_z);
319 
320  _clusters->AddCluster(clusters[cl]);
321 
322  if (Verbosity() > 1)
323  {
324  std::cout << "RawClusterBuilderTopo::export_clusters: added cluster with E = " << clusters_E[cl] << ", eta = " << -1 * log(tan(std::atan2(std::sqrt(mean_y * mean_y + mean_x * mean_x), mean_z) / 2.0)) << ", phi = " << std::atan2(mean_y, mean_x) << std::endl;
325  }
326  }
327 
328  return;
329 }
330 
331 RawClusterBuilderTopo::RawClusterBuilderTopo(const std::string &name)
332  : SubsysReco(name)
333 {
334  // geometry defined at run-time
335  _EMCAL_NETA = -1;
336  _EMCAL_NPHI = -1;
337 
338  _HCAL_NETA = -1;
339  _HCAL_NPHI = -1;
340  std::fill(std::begin(_geom_containers), std::end(_geom_containers), nullptr);
341  _noise_LAYER[0] = 0.0025;
342  _noise_LAYER[1] = 0.006;
343  _noise_LAYER[2] = 0.03; // EM
344 
345  _sigma_seed = 4.0;
346  _sigma_grow = 2.0;
347  _sigma_peri = 0.0;
348 
349  _allow_corner_neighbor = true;
350 
351  _enable_HCal = true;
352  _enable_EMCal = true;
353 
354  _do_split = true;
355  _R_shower = 0.025;
356 
357  _local_max_minE_LAYER[0] = 1;
358  _local_max_minE_LAYER[1] = 1;
359  _local_max_minE_LAYER[2] = 1;
360 
361  ClusterNodeName = "TOPOCLUSTER_HCAL";
362 }
363 
364 int RawClusterBuilderTopo::InitRun(PHCompositeNode *topNode)
365 {
366  try
367  {
368  CreateNodes(topNode);
369  }
370  catch (std::exception &e)
371  {
372  std::cout << PHWHERE << ": " << e.what() << std::endl;
373  throw;
374  }
375 
376  if (Verbosity() > 0)
377  {
378  std::cout << "RawClusterBuilderTopo::InitRun: initialized with EMCal enable = " << _enable_EMCal << " and I+OHCal enable = " << _enable_HCal << std::endl;
379  std::cout << "RawClusterBuilderTopo::InitRun: initialized with sigma_noise in EMCal / IHCal / OHCal = " << _noise_LAYER[2] << " / " << _noise_LAYER[0] << " / " << _noise_LAYER[1] << std::endl;
380  std::cout << "RawClusterBuilderTopo::InitRun: initialized with noise multiples for seeding / growth / perimeter ( S / N / P ) = " << _sigma_seed << " / " << _sigma_grow << " / " << _sigma_peri << std::endl;
381  std::cout << "RawClusterBuilderTopo::InitRun: initialized with allow_corner_neighbor = " << _allow_corner_neighbor << " (in HCal)" << std::endl;
382  std::cout << "RawClusterBuilderTopo::InitRun: initialized with do_split = " << _do_split << " , R_shower = " << _R_shower << " (angular units) " << std::endl;
383  std::cout << "RawClusterBuilderTopo::InitRun: initialized with minE for local max in EMCal / IHCal / OHCal = " << _local_max_minE_LAYER[2] << " / " << _local_max_minE_LAYER[0] << " / " << _local_max_minE_LAYER[1] << std::endl;
384  }
385 
386  return Fun4AllReturnCodes::EVENT_OK;
387 }
388 
389 int RawClusterBuilderTopo::process_event(PHCompositeNode *topNode)
390 {
391  TowerInfoContainer *towerinfosEM = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_CEMC");
392  TowerInfoContainer *towerinfosIH = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_HCALIN");
393  TowerInfoContainer *towerinfosOH = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_HCALOUT");
394 
395  if (!towerinfosEM)
396  {
397  std::cout << " RawClusterBuilderTopo::process_event : container TOWERINFO_CALIB_CEMC does not exist, aborting " << std::endl;
398  return Fun4AllReturnCodes::ABORTEVENT;
399  }
400  if (!towerinfosIH)
401  {
402  std::cout << " RawClusterBuilderTopo::process_event : container TOWERINFO_CALIB_HCALIN does not exist, aborting " << std::endl;
403  return Fun4AllReturnCodes::ABORTEVENT;
404  }
405  if (!towerinfosOH)
406  {
407  std::cout << " RawClusterBuilderTopo::process_event : container TOWERINFO_CALIB_HCALOUT does not exist, aborting " << std::endl;
408  return Fun4AllReturnCodes::ABORTEVENT;
409  }
410 
411  _geom_containers[0] = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALIN");
412  _geom_containers[1] = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALOUT");
413  _geom_containers[2] = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_CEMC");
414 
415  if (!_geom_containers[0])
416  {
417  std::cout << " RawClusterBuilderTopo::process_event : container TOWERGEOM_HCALIN does not exist, aborting " << std::endl;
418  return Fun4AllReturnCodes::ABORTEVENT;
419  }
420  if (!_geom_containers[1])
421  {
422  std::cout << " RawClusterBuilderTopo::process_event : container TOWERGEOM_HCALOUT does not exist, aborting " << std::endl;
423  return Fun4AllReturnCodes::ABORTEVENT;
424  }
425  if (!_geom_containers[2])
426  {
427  std::cout << " RawClusterBuilderTopo::process_event : container TOWERGEOM_CEMC does not exist, aborting " << std::endl;
428  return Fun4AllReturnCodes::ABORTEVENT;
429  }
430 
431  if (Verbosity() > 10)
432  {
433  std::cout << "RawClusterBuilderTopo::process_event: " << towerinfosEM->size() << " TOWERINFO_CALIB_CEMC towers" << std::endl;
434  std::cout << "RawClusterBuilderTopo::process_event: " << towerinfosIH->size() << " TOWERINFO_CALIB_HCALIN towers" << std::endl;
435  std::cout << "RawClusterBuilderTopo::process_event: " << towerinfosOH->size() << " TOWERINFO_CALIB_HCALOUT towers" << std::endl;
436 
437  std::cout << "RawClusterBuilderTopo::process_event: pointer to TOWERGEOM_CEMC: " << _geom_containers[2] << std::endl;
438  std::cout << "RawClusterBuilderTopo::process_event: pointer to TOWERGEOM_HCALIN: " << _geom_containers[0] << std::endl;
439  std::cout << "RawClusterBuilderTopo::process_event: pointer to TOWERGEOM_HCALOUT: " << _geom_containers[1] << std::endl;
440  }
441 
442  if (_EMCAL_NETA < 0)
443  {
444  // define geometry only once if it has not been yet
445  _EMCAL_NETA = _geom_containers[2]->get_etabins();
446  _EMCAL_NPHI = _geom_containers[2]->get_phibins();
447 
448  _EMTOWERMAP_STATUS_ETA_PHI.resize(_EMCAL_NETA, std::vector<int>(_EMCAL_NPHI, -2));
449  _EMTOWERMAP_KEY_ETA_PHI.resize(_EMCAL_NETA, std::vector<int>(_EMCAL_NPHI, 0));
450  _EMTOWERMAP_E_ETA_PHI.resize(_EMCAL_NETA, std::vector<float>(_EMCAL_NPHI, 0));
451  }
452 
453  if (_HCAL_NETA < 0)
454  {
455  // define geometry only once if it has not been yet
456  _HCAL_NETA = _geom_containers[1]->get_etabins();
457  _HCAL_NPHI = _geom_containers[1]->get_phibins();
458 
459  _TOWERMAP_STATUS_LAYER_ETA_PHI.resize(2, std::vector<std::vector<int> >(_HCAL_NETA, std::vector<int>(_HCAL_NPHI, -2)));
460  _TOWERMAP_KEY_LAYER_ETA_PHI.resize(2, std::vector<std::vector<int> >(_HCAL_NETA, std::vector<int>(_HCAL_NPHI, 0)));
461  _TOWERMAP_E_LAYER_ETA_PHI.resize(2, std::vector<std::vector<float> >(_HCAL_NETA, std::vector<float>(_HCAL_NPHI, 0)));
462  }
463 
464  // reset maps
465  // but note -- do not reset keys!
466  for (int ieta = 0; ieta < _EMCAL_NETA; ieta++)
467  {
468  for (int iphi = 0; iphi < _EMCAL_NPHI; iphi++)
469  {
470  _EMTOWERMAP_STATUS_ETA_PHI[ieta][iphi] = -2; // set tower does not exist
471  _EMTOWERMAP_E_ETA_PHI[ieta][iphi] = 0; // set zero energy
472  }
473  }
474  for (int ilayer = 0; ilayer < 2; ilayer++)
475  {
476  for (int ieta = 0; ieta < _HCAL_NETA; ieta++)
477  {
478  for (int iphi = 0; iphi < _HCAL_NPHI; iphi++)
479  {
480  _TOWERMAP_STATUS_LAYER_ETA_PHI[ilayer][ieta][iphi] = -2; // set tower does not exist
481  _TOWERMAP_E_LAYER_ETA_PHI[ilayer][ieta][iphi] = 0; // set zero energy
482  }
483  }
484  }
485 
486  // setup
487  std::vector<std::pair<int, float> > list_of_seeds;
488 
489  // translate towers to our internal representation
490  if (_enable_EMCal)
491  {
492  TowerInfo *towerInfo = nullptr;
493  for (unsigned int iEM = 0; iEM < towerinfosEM->size(); iEM++)
494  {
495  towerInfo = towerinfosEM->get_tower_at_channel(iEM);
496  unsigned int towerinfo_key = towerinfosEM->encode_key(iEM);
497  int ti_ieta = towerinfosEM->getTowerEtaBin(towerinfo_key);
498  int ti_iphi = towerinfosEM->getTowerPhiBin(towerinfo_key);
499  const RawTowerDefs::keytype key = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::CEMC, ti_ieta, ti_iphi);
500 
501  RawTowerGeom *tower_geom = _geom_containers[2]->get_tower_geometry(key);
502 
503  int ieta = _geom_containers[2]->get_etabin(tower_geom->get_eta());
504  int iphi = _geom_containers[2]->get_phibin(tower_geom->get_phi());
505  float this_E = towerInfo->get_energy();
506 
507  if (this_E < 1.E-10)
508  {
509  continue;
510  }
511 
512  _EMTOWERMAP_STATUS_ETA_PHI[ieta][iphi] = -1; // change status to unknown
513  _EMTOWERMAP_E_ETA_PHI[ieta][iphi] = this_E;
514  _EMTOWERMAP_KEY_ETA_PHI[ieta][iphi] = key;
515 
516  if (this_E > _sigma_seed * _noise_LAYER[2])
517  {
518  int ID = get_ID(2, ieta, iphi);
519  list_of_seeds.emplace_back(ID, this_E);
520  if (Verbosity() > 10)
521  {
522  std::cout << "RawClusterBuilderTopo::process_event: adding EMCal tower at ieta / iphi = " << ieta << " / " << iphi << " with E = " << this_E << std::endl;
523  std::cout << " --> ID = " << ID << " , check ilayer / ieta / iphi = " << get_ilayer_from_ID(ID) << " / " << get_ieta_from_ID(ID) << " / " << get_iphi_from_ID(ID) << std::endl;
524  };
525  }
526  }
527  }
528 
529  // translate towers to our internal representation
530  if (_enable_HCal)
531  {
532  TowerInfo *towerInfo = nullptr;
533  for (unsigned int iIH = 0; iIH < towerinfosIH->size(); iIH++)
534  {
535  towerInfo = towerinfosIH->get_tower_at_channel(iIH);
536  unsigned int towerinfo_key = towerinfosIH->encode_key(iIH);
537  int ti_ieta = towerinfosIH->getTowerEtaBin(towerinfo_key);
538  int ti_iphi = towerinfosIH->getTowerPhiBin(towerinfo_key);
539  const RawTowerDefs::keytype key = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::HCALIN, ti_ieta, ti_iphi);
540 
541  RawTowerGeom *tower_geom = _geom_containers[0]->get_tower_geometry(key);
542 
543  int ieta = _geom_containers[0]->get_etabin(tower_geom->get_eta());
544  int iphi = _geom_containers[0]->get_phibin(tower_geom->get_phi());
545  float this_E = towerInfo->get_energy();
546 
547  if (this_E < 1.E-10)
548  {
549  continue;
550  }
551 
552  _TOWERMAP_STATUS_LAYER_ETA_PHI[0][ieta][iphi] = -1; // change status to unknown
553  _TOWERMAP_E_LAYER_ETA_PHI[0][ieta][iphi] = this_E;
554  _TOWERMAP_KEY_LAYER_ETA_PHI[0][ieta][iphi] = key;
555 
556  if (this_E > _sigma_seed * _noise_LAYER[0])
557  {
558  int ID = get_ID(0, ieta, iphi);
559  list_of_seeds.emplace_back(ID, this_E);
560  if (Verbosity() > 10)
561  {
562  std::cout << "RawClusterBuilderTopo::process_event: adding IHCal tower at ieta / iphi = " << ieta << " / " << iphi << " with E = " << this_E << std::endl;
563  std::cout << " --> ID = " << ID << " , check ilayer / ieta / iphi = " << get_ilayer_from_ID(ID) << " / " << get_ieta_from_ID(ID) << " / " << get_iphi_from_ID(ID) << std::endl;
564  };
565  }
566  }
567 
568  for (unsigned int iOH = 0; iOH < towerinfosOH->size(); iOH++)
569  {
570  towerInfo = towerinfosOH->get_tower_at_channel(iOH);
571  unsigned int towerinfo_key = towerinfosOH->encode_key(iOH);
572  int ti_ieta = towerinfosOH->getTowerEtaBin(towerinfo_key);
573  int ti_iphi = towerinfosOH->getTowerPhiBin(towerinfo_key);
574  const RawTowerDefs::keytype key = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::HCALOUT, ti_ieta, ti_iphi);
575 
576  RawTowerGeom *tower_geom = _geom_containers[1]->get_tower_geometry(key);
577 
578  int ieta = _geom_containers[1]->get_etabin(tower_geom->get_eta());
579  int iphi = _geom_containers[1]->get_phibin(tower_geom->get_phi());
580  float this_E = towerInfo->get_energy();
581 
582  if (this_E < 1.E-10)
583  {
584  continue;
585  }
586 
587  _TOWERMAP_STATUS_LAYER_ETA_PHI[1][ieta][iphi] = -1; // change status to unknown
588  _TOWERMAP_E_LAYER_ETA_PHI[1][ieta][iphi] = this_E;
589  _TOWERMAP_KEY_LAYER_ETA_PHI[1][ieta][iphi] = key;
590 
591  if (this_E > _sigma_seed * _noise_LAYER[1])
592  {
593  int ID = get_ID(1, ieta, iphi);
594  list_of_seeds.emplace_back(ID, this_E);
595  if (Verbosity() > 10)
596  {
597  std::cout << "RawClusterBuilderTopo::process_event: adding OHCal tower at ieta / iphi = " << ieta << " / " << iphi << " with E = " << this_E << std::endl;
598  std::cout << " --> ID = " << ID << " , check ilayer / ieta / iphi = " << get_ilayer_from_ID(ID) << " / " << get_ieta_from_ID(ID) << " / " << get_iphi_from_ID(ID) << std::endl;
599  };
600  }
601  }
602  }
603 
604  if (Verbosity() > 10)
605  {
606  for (unsigned int n = 0; n < list_of_seeds.size(); n++)
607  {
608  std::cout << "RawClusterBuilderTopo::process_event: unsorted seed element n = " << n << " , ID / E = " << list_of_seeds.at(n).first << " / " << list_of_seeds.at(n).second << std::endl;
609  }
610  }
611 
612  std::sort(list_of_seeds.begin(), list_of_seeds.end(), sort_by_pair_second);
613 
614  if (Verbosity() > 10)
615  {
616  for (unsigned int n = 0; n < list_of_seeds.size(); n++)
617  {
618  std::cout << "RawClusterBuilderTopo::process_event: sorted seed element n = " << n << " , ID / E = " << list_of_seeds.at(n).first << " / " << list_of_seeds.at(n).second << std::endl;
619  }
620  }
621 
622  if (Verbosity() > 0)
623  {
624  std::cout << "RawClusterBuilderTopo::process_event: initialized with " << list_of_seeds.size() << " seeds with E > 4*sigma " << std::endl;
625  }
626 
627  int cluster_index = 0; // begin counting clusters
628 
629  std::vector<std::vector<int> > all_cluster_towers; // store final cluster tower lists here
630 
631  while (list_of_seeds.size() > 0)
632  {
633  int seed_ID = list_of_seeds.at(0).first;
634  list_of_seeds.erase(list_of_seeds.begin());
635 
636  if (Verbosity() > 5)
637  {
638  std::cout << " RawClusterBuilderTopo::process_event: in seeded loop, current seed has ID = " << seed_ID << " , length of remaining seed vector = " << list_of_seeds.size() << std::endl;
639  }
640 
641  // if this seed was already claimed by some other seed during its growth, remove it and do nothing
642  int seed_status = get_status_from_ID(seed_ID);
643  if (seed_status > -1)
644  {
645  if (Verbosity() > 10)
646  {
647  std::cout << " --> already owned by cluster # " << seed_status << std::endl;
648  }
649  continue; // go onto the next iteration of the loop
650  }
651 
652  // this seed tower now owned by new cluster
653  set_status_by_ID(seed_ID, cluster_index);
654 
655  std::vector<int> cluster_tower_ID;
656  cluster_tower_ID.push_back(seed_ID);
657 
658  std::vector<int> grow_tower_ID;
659  grow_tower_ID.push_back(seed_ID);
660 
661  // iteratively process growth towers, adding > 2 * sigma neighbors to the list for further checking
662 
663  if (Verbosity() > 5)
664  {
665  std::cout << " RawClusterBuilderTopo::process_event: Entering Growth stage for cluster " << cluster_index << std::endl;
666  }
667 
668  while (grow_tower_ID.size() > 0)
669  {
670  int grow_ID = grow_tower_ID.at(0);
671  grow_tower_ID.erase(grow_tower_ID.begin());
672 
673  if (Verbosity() > 5)
674  {
675  std::cout << " --> cluster " << cluster_index << ", growth stage, examining neighbors of ID " << grow_ID << ", " << grow_tower_ID.size() << " grow towers left" << std::endl;
676  }
677 
678  std::vector<int> adjacent_tower_IDs = get_adjacent_towers_by_ID(grow_ID);
679 
680  for (int this_adjacent_tower_ID : adjacent_tower_IDs)
681  {
682  if (Verbosity() > 10)
683  {
684  std::cout << " --> --> --> checking possible adjacent tower with ID " << this_adjacent_tower_ID << " : ";
685  }
686  int test_layer = get_ilayer_from_ID(this_adjacent_tower_ID);
687 
688  // if tower does not exist, continue
689  if (get_status_from_ID(this_adjacent_tower_ID) == -2)
690  {
691  if (Verbosity() > 10)
692  {
693  std::cout << "does not exist " << std::endl;
694  }
695  continue;
696  }
697 
698  // if tower is owned by THIS cluster already, continue
699  if (get_status_from_ID(this_adjacent_tower_ID) == cluster_index)
700  {
701  if (Verbosity() > 10)
702  {
703  std::cout << "already owned by this cluster index " << cluster_index << std::endl;
704  }
705  continue;
706  }
707 
708  // if tower has < 2*sigma energy, continue
709  if (get_E_from_ID(this_adjacent_tower_ID) < _sigma_grow * _noise_LAYER[test_layer])
710  {
711  if (Verbosity() > 10)
712  {
713  std::cout << "E = " << get_E_from_ID(this_adjacent_tower_ID) << " under 2*sigma threshold " << std::endl;
714  }
715  continue;
716  }
717 
718  // if tower is owned by somebody else, continue (although should this really happen?)
719  if (get_status_from_ID(this_adjacent_tower_ID) > -1)
720  {
721  if (Verbosity() > 10)
722  {
723  std::cout << "ERROR! in growth stage, encountered >2sigma tower which is already owned?!" << std::endl;
724  }
725  continue;
726  }
727 
728  // tower good to be added to cluster and to list of grow towers
729  grow_tower_ID.push_back(this_adjacent_tower_ID);
730  cluster_tower_ID.push_back(this_adjacent_tower_ID);
731  set_status_by_ID(this_adjacent_tower_ID, cluster_index);
732  if (Verbosity() > 10)
733  {
734  std::cout << "add this tower ( ID " << this_adjacent_tower_ID << " ) to grow list " << std::endl;
735  }
736  }
737 
738  if (Verbosity() > 5)
739  {
740  std::cout << " --> after examining neighbors, grow list is now " << grow_tower_ID.size() << ", # of towers in cluster = " << cluster_tower_ID.size() << std::endl;
741  }
742  }
743 
744  // done growing cluster, now add on perimeter towers with E > 0 * sigma
745  if (Verbosity() > 5)
746  {
747  std::cout << " RawClusterBuilderTopo::process_event: Entering Perimeter stage for cluster " << cluster_index << std::endl;
748  }
749  // we'll be adding on to the cluster list, so get the # of core towers first
750  int n_core_towers = cluster_tower_ID.size();
751 
752  for (int ic = 0; ic < n_core_towers; ic++)
753  {
754  int core_ID = cluster_tower_ID.at(ic);
755 
756  if (Verbosity() > 5)
757  {
758  std::cout << " --> cluster " << cluster_index << ", perimeter stage, examining neighbors of ID " << core_ID << ", core cluster # " << ic << " of " << n_core_towers << " total " << std::endl;
759  }
760  std::vector<int> adjacent_tower_IDs = get_adjacent_towers_by_ID(core_ID);
761 
762  for (int this_adjacent_tower_ID : adjacent_tower_IDs)
763  {
764  if (Verbosity() > 10)
765  {
766  std::cout << " --> --> --> checking possible adjacent tower with ID " << this_adjacent_tower_ID << " : ";
767  }
768 
769  int test_layer = get_ilayer_from_ID(this_adjacent_tower_ID);
770 
771  // if tower does not exist, continue
772  if (get_status_from_ID(this_adjacent_tower_ID) == -2)
773  {
774  if (Verbosity() > 10)
775  {
776  std::cout << "does not exist " << std::endl;
777  }
778  continue;
779  }
780 
781  // if tower is owned by somebody else (including current cluster), continue. ( allowed during perimeter fixing state )
782  if (get_status_from_ID(this_adjacent_tower_ID) > -1)
783  {
784  if (Verbosity() > 10)
785  {
786  std::cout << "already owned by other cluster index " << get_status_from_ID(this_adjacent_tower_ID) << std::endl;
787  }
788  continue;
789  }
790 
791  // if tower has < 0*sigma energy, continue
792  if (get_E_from_ID(this_adjacent_tower_ID) < _sigma_peri * _noise_LAYER[test_layer])
793  {
794  if (Verbosity() > 10)
795  {
796  std::cout << "E = " << get_E_from_ID(this_adjacent_tower_ID) << " under 0*sigma threshold " << std::endl;
797  }
798  continue;
799  }
800 
801  // perimeter tower good to be added to cluster
802  cluster_tower_ID.push_back(this_adjacent_tower_ID);
803  set_status_by_ID(this_adjacent_tower_ID, cluster_index);
804  if (Verbosity() > 10)
805  {
806  std::cout << "add this tower ( ID " << this_adjacent_tower_ID << " ) to cluster " << std::endl;
807  }
808  }
809 
810  if (Verbosity() > 5)
811  {
812  std::cout << " --> after examining perimeter neighbors, # of towers in cluster is now = " << cluster_tower_ID.size() << std::endl;
813  }
814  }
815 
816  // keep track of these
817  all_cluster_towers.push_back(cluster_tower_ID);
818 
819  // increment cluster index for next one
820  cluster_index++;
821  }
822 
823  if (Verbosity() > 0)
824  {
825  std::cout << "RawClusterBuilderTopo::process_event: " << cluster_index << " topo-clusters initially reconstructed, entering splitting step" << std::endl;
826  }
827 
828  int original_cluster_index = cluster_index; // since it may be updated
829 
830  // now entering cluster splitting stage
831 
832  for (int cl = 0; cl < original_cluster_index; cl++)
833  {
834  std::vector<int> original_towers = all_cluster_towers.at(cl);
835 
836  if (!_do_split)
837  {
838  // don't run splitting, just export entire cluster as it is
839  if (Verbosity() > 2)
840  {
841  std::cout << "RawClusterBuilderTopo::process_event: splitting step disabled, cluster " << cluster_index << " is final" << std::endl;
842  }
843  export_single_cluster(original_towers);
844  continue;
845  }
846 
847  std::vector<std::pair<int, float> > local_maxima_ID;
848 
849  // iterate through each tower, looking for maxima
850  for (int tower_ID : original_towers)
851  {
852  if (Verbosity() > 10)
853  {
854  std::cout << " -> examining tower ID " << tower_ID << " for possible local maximum " << std::endl;
855  }
856 
857  // check minimum energy
858  if (get_E_from_ID(tower_ID) < _local_max_minE_LAYER[get_ilayer_from_ID(tower_ID)])
859  {
860  if (Verbosity() > 10)
861  {
862  std::cout << " -> -> energy E = " << get_E_from_ID(tower_ID) << " < " << _local_max_minE_LAYER[get_ilayer_from_ID(tower_ID)] << " too low" << std::endl;
863  }
864  continue;
865  }
866 
867  // examine neighbors
868  std::vector<int> adjacent_tower_IDs = get_adjacent_towers_by_ID(tower_ID);
869  int neighbors_in_cluster = 0;
870 
871  // check for higher neighbox
872  bool has_higher_neighbor = false;
873  for (int this_adjacent_tower_ID : adjacent_tower_IDs)
874  {
875  if (get_status_from_ID(this_adjacent_tower_ID) != cl)
876  {
877  continue; // only consider neighbors in cluster, obviously
878  }
879 
880  neighbors_in_cluster++;
881 
882  if (get_E_from_ID(this_adjacent_tower_ID) > get_E_from_ID(tower_ID))
883  {
884  if (Verbosity() > 10)
885  {
886  std::cout << " -> -> has higher-energy neighbor ID / E = " << this_adjacent_tower_ID << " / " << get_E_from_ID(this_adjacent_tower_ID) << std::endl;
887  }
888  has_higher_neighbor = true; // at this point we can break -- we won't need to count the number of good neighbors, since we won't even pass the E_neighbor test
889  break;
890  }
891  }
892 
893  if (has_higher_neighbor)
894  {
895  continue; // if we broke out, now continue
896  }
897 
898  // check number of neighbors
899  if (neighbors_in_cluster < 4)
900  {
901  if (Verbosity() > 10)
902  {
903  std::cout << " -> -> too few neighbors N = " << neighbors_in_cluster << std::endl;
904  }
905  continue;
906  }
907 
908  local_maxima_ID.emplace_back(tower_ID, get_E_from_ID(tower_ID));
909  }
910 
911  // check for possible EMCal-OHCal seed overlaps
912  for (unsigned int n = 0; n < local_maxima_ID.size(); n++)
913  {
914  // only look at I/OHCal local maxima
915  std::pair<int, float> this_LM = local_maxima_ID.at(n);
916  if (get_ilayer_from_ID(this_LM.first) == 2)
917  {
918  continue;
919  }
920 
921  float this_phi = _geom_containers[get_ilayer_from_ID(this_LM.first)]->get_phicenter(get_iphi_from_ID(this_LM.first));
922  if (this_phi > M_PI)
923  {
924  this_phi -= 2 * M_PI;
925  }
926  float this_eta = _geom_containers[get_ilayer_from_ID(this_LM.first)]->get_etacenter(get_ieta_from_ID(this_LM.first));
927 
928  bool has_EM_overlap = false;
929 
930  // check all other local maxima for overlaps
931  for (unsigned int n2 = 0; n2 < local_maxima_ID.size(); n2++)
932  {
933  if (n == n2)
934  {
935  continue; // don't check the same one
936  }
937 
938  // only look at EMCal local mazima
939  std::pair<int, float> this_LM2 = local_maxima_ID.at(n2);
940  if (get_ilayer_from_ID(this_LM2.first) != 2)
941  {
942  continue;
943  }
944 
945  float this_phi2 = _geom_containers[get_ilayer_from_ID(this_LM2.first)]->get_phicenter(get_iphi_from_ID(this_LM2.first));
946  if (this_phi2 > M_PI)
947  {
948  this_phi -= 2 * M_PI;
949  }
950  float this_eta2 = _geom_containers[get_ilayer_from_ID(this_LM2.first)]->get_etacenter(get_ieta_from_ID(this_LM2.first));
951 
952  // calculate geometric dR
953  float dR = calculate_dR(this_eta, this_eta2, this_phi, this_phi2);
954 
955  // check for and report overlaps
956  if (dR < 0.15)
957  {
958  has_EM_overlap = true;
959  if (Verbosity() > 2)
960  {
961  std::cout << "RawClusterBuilderTopo::process_event : removing I/OHal local maximum (ID,E,phi,eta = " << this_LM.first << ", " << this_LM.second << ", " << this_phi << ", " << this_eta << "), ";
962  std::cout << "due to EM overlap (ID,E,phi,eta = " << this_LM2.first << ", " << this_LM2.second << ", " << this_phi2 << ", " << this_eta2 << "), dR = " << dR << std::endl;
963  }
964  break;
965  }
966  }
967 
968  if (has_EM_overlap)
969  {
970  // remove the I/OHCal local maximum from the list
971  local_maxima_ID.erase(local_maxima_ID.begin() + n);
972  // make sure to back up one index...
973  n = n - 1;
974  } // otherwise, keep this local maximum
975  }
976 
977  // only now print out full set of local maxima
978  if (Verbosity() > 2)
979  {
980  for (auto this_LM : local_maxima_ID)
981  {
982  int tower_ID = this_LM.first;
983  std::cout << "RawClusterBuilderTopo::process_event in cluster " << cl << ", tower ID " << tower_ID << " is LOCAL MAXIMUM with layer / E = " << get_ilayer_from_ID(tower_ID) << " / " << get_E_from_ID(tower_ID) << ", ";
984  float this_phi = _geom_containers[get_ilayer_from_ID(tower_ID)]->get_phicenter(get_iphi_from_ID(tower_ID));
985  if (this_phi > M_PI)
986  {
987  this_phi -= 2 * M_PI;
988  }
989  std::cout << " eta / phi = " << _geom_containers[get_ilayer_from_ID(tower_ID)]->get_etacenter(get_ieta_from_ID(tower_ID)) << " / " << this_phi << std::endl;
990  }
991  }
992 
993  // do we have only 1 or 0 local maxima?
994  if (local_maxima_ID.size() <= 1)
995  {
996  if (Verbosity() > 2)
997  {
998  std::cout << "RawClusterBuilderTopo::process_event cluster " << cl << " has only " << local_maxima_ID.size() << " local maxima, not splitting " << std::endl;
999  }
1000  export_single_cluster(original_towers);
1001 
1002  continue;
1003  }
1004 
1005  // engage splitting procedure!
1006 
1007  if (Verbosity() > 2)
1008  {
1009  std::cout << "RawClusterBuilderTopo::process_event splitting cluster " << cl << " into " << local_maxima_ID.size() << " according to local maxima!" << std::endl;
1010  }
1011  // translate all cluster towers to a map which keeps track of their ownership
1012  // -1 means unseen
1013  // -2 means seen and in the seed list now (e.g. don't add it to the seed list again)
1014  // -3 shared tower, ignore going forward...
1015  std::map<int, std::pair<int, int> > tower_ownership;
1016  for (int &original_tower : original_towers)
1017  {
1018  tower_ownership[original_tower] = std::pair<int, int>(-1, -1); // initialize all towers as un-seen
1019  }
1020  std::vector<int> seed_list;
1021  std::vector<int> neighbor_list;
1022  std::vector<int> shared_list;
1023 
1024  // sort maxima before populating seed list
1025  std::sort(local_maxima_ID.begin(), local_maxima_ID.end(), sort_by_pair_second);
1026 
1027  // initialize neighbor list
1028  for (unsigned int s = 0; s < local_maxima_ID.size(); s++)
1029  {
1030  tower_ownership[local_maxima_ID.at(s).first] = std::pair<int, int>(s, -1);
1031  neighbor_list.push_back(local_maxima_ID.at(s).first);
1032  }
1033 
1034  if (Verbosity() > 100)
1035  {
1036  for (int &original_tower : original_towers)
1037  {
1038  std::pair<int, int> the_pair = tower_ownership[original_tower];
1039  std::cout << " Debug Pre-Split: tower_ownership[ " << original_tower << " ] = ( " << the_pair.first << ", " << the_pair.second << " ) ";
1040  std::cout << " , layer / ieta / iphi = " << get_ilayer_from_ID(original_tower) << " / " << get_ieta_from_ID(original_tower) << " / " << get_iphi_from_ID(original_tower);
1041  std::cout << std::endl;
1042  }
1043  }
1044 
1045  bool first_pass = true;
1046 
1047  do
1048  {
1049  if (Verbosity() > 5)
1050  {
1051  std::cout << " -> starting split loop with " << seed_list.size() << " seed, " << neighbor_list.size() << " neighbor, and " << shared_list.size() << " shared towers " << std::endl;
1052  }
1053  // go through neighbor list, assigning ownership only via the seed list
1054  std::vector<int> new_ownerships;
1055 
1056  for (unsigned int n = 0; n < neighbor_list.size(); n++)
1057  {
1058  int neighbor_ID = neighbor_list.at(n);
1059 
1060  if (Verbosity() > 10)
1061  {
1062  std::cout << " -> -> looking at neighbor " << n << " (tower ID " << neighbor_ID << " ) of " << neighbor_list.size() << " total" << std::endl;
1063  }
1064  if (first_pass)
1065  {
1066  if (Verbosity() > 10)
1067  {
1068  std::cout << " -> -> -> special first pass rules, this tower already owned by pseudocluster " << tower_ownership[neighbor_ID].first << std::endl;
1069  }
1070  new_ownerships.push_back(tower_ownership[neighbor_ID].first);
1071  }
1072  else
1073  {
1074  std::map<int, bool> pseudocluster_adjacency;
1075  for (unsigned int s = 0; s < local_maxima_ID.size(); s++)
1076  {
1077  pseudocluster_adjacency[s] = false;
1078  }
1079  // look over all towers THIS one is adjacent to, and count up...
1080  std::vector<int> adjacent_tower_IDs = get_adjacent_towers_by_ID(neighbor_ID);
1081 
1082  for (int this_adjacent_tower_ID : adjacent_tower_IDs)
1083  {
1084  if (get_status_from_ID(this_adjacent_tower_ID) != cl)
1085  {
1086  continue;
1087  }
1088 
1089  if (tower_ownership[this_adjacent_tower_ID].first > -1)
1090  {
1091  if (Verbosity() > 20)
1092  {
1093  std::cout << " -> -> -> adjacent tower to this one, with ID " << this_adjacent_tower_ID << " , is owned by pseudocluster " << tower_ownership[this_adjacent_tower_ID].first << std::endl;
1094  }
1095  pseudocluster_adjacency[tower_ownership[this_adjacent_tower_ID].first] = true;
1096  }
1097  }
1098  int n_pseudocluster_adjacent = 0;
1099  int last_adjacent_pseudocluster = -1;
1100  for (unsigned int s = 0; s < local_maxima_ID.size(); s++)
1101  {
1102  if (pseudocluster_adjacency[s])
1103  {
1104  last_adjacent_pseudocluster = s;
1105  n_pseudocluster_adjacent++;
1106  if (Verbosity() > 20)
1107  {
1108  std::cout << " -> -> adjacent to pseudocluster " << s << std::endl;
1109  }
1110  }
1111  }
1112 
1113  if (n_pseudocluster_adjacent == 0)
1114  {
1115  std::cout << " -> -> ERROR! How can a neighbor tower at this stage be adjacent to no pseudoclusters?? " << std::endl;
1116  new_ownerships.push_back(9999);
1117  }
1118  else if (n_pseudocluster_adjacent == 1)
1119  {
1120  if (Verbosity() > 10)
1121  {
1122  std::cout << " -> -> neighbor tower " << neighbor_ID << " is ONLY adjacent to one pseudocluster # " << last_adjacent_pseudocluster << std::endl;
1123  }
1124  new_ownerships.push_back(last_adjacent_pseudocluster);
1125  }
1126  else
1127  {
1128  if (Verbosity() > 10)
1129  {
1130  std::cout << " -> -> neighbor tower " << neighbor_ID << " is adjacent to " << n_pseudocluster_adjacent << " pseudoclusters, move to shared list " << std::endl;
1131  }
1132  new_ownerships.push_back(-3);
1133  }
1134  }
1135  }
1136 
1137  if (Verbosity() > 5)
1138  {
1139  std::cout << " -> now updating status of all " << neighbor_list.size() << " original neighbors " << std::endl;
1140  }
1141  // transfer neighbor list to seed list or shared list
1142  for (unsigned int n = 0; n < neighbor_list.size(); n++)
1143  {
1144  int neighbor_ID = neighbor_list.at(n);
1145  if (new_ownerships.at(n) > -1)
1146  {
1147  tower_ownership[neighbor_ID] = std::pair<int, int>(new_ownerships.at(n), -1);
1148  seed_list.push_back(neighbor_ID);
1149  if (Verbosity() > 20)
1150  {
1151  std::cout << " -> -> neighbor ID " << neighbor_ID << " has new status " << new_ownerships.at(n) << std::endl;
1152  }
1153  }
1154  if (new_ownerships.at(n) == -3)
1155  {
1156  tower_ownership[neighbor_ID] = std::pair<int, int>(-3, -1);
1157  shared_list.push_back(neighbor_ID);
1158  if (Verbosity() > 20)
1159  {
1160  std::cout << " -> -> neighbor ID " << neighbor_ID << " has new status " << -3 << std::endl;
1161  }
1162  }
1163  }
1164 
1165  if (Verbosity() > 5)
1166  {
1167  std::cout << " producing a new neighbor list ... " << std::endl;
1168  }
1169  // populate a new neighbor list from the about-to-be-owned towers before transferring this one
1170  std::list<int> new_neighbor_list;
1171  for (unsigned int n = 0; n < neighbor_list.size(); n++)
1172  {
1173  int neighbor_ID = neighbor_list.at(n);
1174  if (new_ownerships.at(n) > -1)
1175  {
1176  std::vector<int> adjacent_tower_IDs = get_adjacent_towers_by_ID(neighbor_ID);
1177 
1178  for (int this_adjacent_tower_ID : adjacent_tower_IDs)
1179  {
1180  if (get_status_from_ID(this_adjacent_tower_ID) != cl)
1181  {
1182  continue;
1183  }
1184  if (tower_ownership[this_adjacent_tower_ID].first == -1)
1185  {
1186  new_neighbor_list.push_back(this_adjacent_tower_ID);
1187  if (Verbosity() > 5)
1188  {
1189  std::cout << " -> queueing up to add tower " << this_adjacent_tower_ID << " , neighbor of tower " << neighbor_ID << " to new neighbor list" << std::endl;
1190  }
1191  }
1192  }
1193  }
1194  }
1195 
1196  if (Verbosity() > 5)
1197  {
1198  std::cout << " new neighbor list has size " << new_neighbor_list.size() << ", but after removing duplicate elements: ";
1199  new_neighbor_list.sort();
1200  new_neighbor_list.unique();
1201  std::cout << new_neighbor_list.size() << std::endl;
1202  }
1203 
1204  // clear neighbor list!
1205  neighbor_list.clear();
1206 
1207  // now transfer over new neighbor list
1208  for (; new_neighbor_list.size() > 0;)
1209  {
1210  neighbor_list.push_back(new_neighbor_list.front());
1211  new_neighbor_list.pop_front();
1212  }
1213 
1214  first_pass = false;
1215 
1216  } while (neighbor_list.size() > 0);
1217 
1218  if (Verbosity() > 100)
1219  {
1220  for (int &original_tower : original_towers)
1221  {
1222  std::pair<int, int> the_pair = tower_ownership[original_tower];
1223  std::cout << " Debug Mid-Split: tower_ownership[ " << original_tower << " ] = ( " << the_pair.first << ", " << the_pair.second << " ) ";
1224  std::cout << " , layer / ieta / iphi = " << get_ilayer_from_ID(original_tower) << " / " << get_ieta_from_ID(original_tower) << " / " << get_iphi_from_ID(original_tower);
1225  std::cout << std::endl;
1226  if (the_pair.first == -1)
1227  {
1228  std::vector<int> adjacent_tower_IDs = get_adjacent_towers_by_ID(original_tower);
1229 
1230  for (int this_adjacent_tower_ID : adjacent_tower_IDs)
1231  {
1232  if (get_status_from_ID(this_adjacent_tower_ID) != cl)
1233  {
1234  continue;
1235  }
1236  std::cout << " -> adjacent to add tower " << this_adjacent_tower_ID << " , which has status " << tower_ownership[this_adjacent_tower_ID].first << std::endl;
1237  }
1238  }
1239  }
1240  }
1241 
1242  // calculate pseudocluster energies and positions
1243  std::vector<float> pseudocluster_sumeta;
1244  std::vector<float> pseudocluster_sumphi;
1245  std::vector<float> pseudocluster_sumE;
1246  std::vector<int> pseudocluster_ntower;
1247  std::vector<float> pseudocluster_eta;
1248  std::vector<float> pseudocluster_phi;
1249 
1250  pseudocluster_sumeta.resize(local_maxima_ID.size(), 0);
1251  pseudocluster_sumphi.resize(local_maxima_ID.size(), 0);
1252  pseudocluster_sumE.resize(local_maxima_ID.size(), 0);
1253  pseudocluster_ntower.resize(local_maxima_ID.size(), 0);
1254 
1255  for (int &original_tower : original_towers)
1256  {
1257  std::pair<int, int> the_pair = tower_ownership[original_tower];
1258  if (the_pair.first > -1)
1259  {
1260  float this_ID = original_tower;
1261  pseudocluster_sumE[the_pair.first] += get_E_from_ID(this_ID);
1262  float this_eta = _geom_containers[get_ilayer_from_ID(this_ID)]->get_etacenter(get_ieta_from_ID(this_ID));
1263  float this_phi = _geom_containers[get_ilayer_from_ID(this_ID)]->get_phicenter(get_iphi_from_ID(this_ID));
1264  // float this_phi = ( get_ilayer_from_ID( this_ID ) == 2 ? geomEM->get_phicenter( get_iphi_from_ID( this_ID ) ) : geomOH->get_phicenter( get_iphi_from_ID( this_ID ) ) );
1265  pseudocluster_sumeta[the_pair.first] += this_eta;
1266  pseudocluster_sumphi[the_pair.first] += this_phi;
1267  pseudocluster_ntower[the_pair.first] += 1;
1268  }
1269  }
1270 
1271  for (unsigned int pc = 0; pc < local_maxima_ID.size(); pc++)
1272  {
1273  pseudocluster_eta.push_back(pseudocluster_sumeta.at(pc) / pseudocluster_ntower.at(pc));
1274  pseudocluster_phi.push_back(pseudocluster_sumphi.at(pc) / pseudocluster_ntower.at(pc));
1275 
1276  if (Verbosity() > 2)
1277  {
1278  std::cout << "RawClusterBuilderTopo::process_event pseudocluster #" << pc << ", E / eta / phi / Ntower = " << pseudocluster_sumE.at(pc) << " / " << pseudocluster_eta.at(pc) << " / " << pseudocluster_phi.at(pc) << " / " << pseudocluster_ntower.at(pc) << std::endl;
1279  }
1280  }
1281 
1282  if (Verbosity() > 2)
1283  {
1284  std::cout << "RawClusterBuilderTopo::process_event now splitting up shared clusters (including unassigned clusters), initial shared list has size " << shared_list.size() << std::endl;
1285  }
1286  // iterate through shared cells, identifying which two they belong to
1287  while (shared_list.size() > 0)
1288  {
1289  // pick the first cell and pop off list
1290  int shared_ID = shared_list.at(0);
1291  shared_list.erase(shared_list.begin());
1292 
1293  if (Verbosity() > 5)
1294  {
1295  std::cout << " -> looking at shared tower " << shared_ID << ", after this one there are " << shared_list.size() << " shared towers left " << std::endl;
1296  }
1297  // look through adjacent pseudoclusters, taking two with highest energies
1298  std::vector<bool> pseudocluster_adjacency;
1299  pseudocluster_adjacency.resize(local_maxima_ID.size(), false);
1300 
1301  std::vector<int> adjacent_tower_IDs = get_adjacent_towers_by_ID(shared_ID);
1302 
1303  for (int this_adjacent_tower_ID : adjacent_tower_IDs)
1304  {
1305  if (get_status_from_ID(this_adjacent_tower_ID) != cl)
1306  {
1307  continue;
1308  }
1309  if (tower_ownership[this_adjacent_tower_ID].first > -1)
1310  {
1311  pseudocluster_adjacency[tower_ownership[this_adjacent_tower_ID].first] = true;
1312  }
1313  if (tower_ownership[this_adjacent_tower_ID].second > -1)
1314  { // can inherit adjacency from shared cluster
1315  pseudocluster_adjacency[tower_ownership[this_adjacent_tower_ID].second] = true;
1316  }
1317  // at the same time, add unowned towers to the list for later examination
1318  if (tower_ownership[this_adjacent_tower_ID].first == -1)
1319  {
1320  shared_list.push_back(this_adjacent_tower_ID);
1321  tower_ownership[this_adjacent_tower_ID] = std::pair<int, int>(-3, -1);
1322  if (Verbosity() > 10)
1323  {
1324  std::cout << " -> while looking at neighbors, have added un-examined tower " << this_adjacent_tower_ID << " to shared list " << std::endl;
1325  }
1326  }
1327  }
1328 
1329  // now figure out which pseudoclustes this shared tower is adjacent to...
1330  int highest_pseudocluster_index = -1;
1331  int second_highest_pseudocluster_index = -1;
1332 
1333  float highest_pseudocluster_E = -1;
1334  float second_highest_pseudocluster_E = -2;
1335 
1336  for (unsigned int n = 0; n < pseudocluster_adjacency.size(); n++)
1337  {
1338  if (!pseudocluster_adjacency[n])
1339  {
1340  continue;
1341  }
1342 
1343  if (pseudocluster_sumE[n] > highest_pseudocluster_E)
1344  {
1345  second_highest_pseudocluster_E = highest_pseudocluster_E;
1346  second_highest_pseudocluster_index = highest_pseudocluster_index;
1347 
1348  highest_pseudocluster_E = pseudocluster_sumE[n];
1349  highest_pseudocluster_index = n;
1350  }
1351  else if (pseudocluster_sumE[n] > second_highest_pseudocluster_E)
1352  {
1353  second_highest_pseudocluster_E = pseudocluster_sumE[n];
1354  second_highest_pseudocluster_index = n;
1355  }
1356  }
1357 
1358  if (Verbosity() > 5)
1359  {
1360  std::cout << " -> highest pseudoclusters its adjacent to are " << highest_pseudocluster_index << " ( E = " << highest_pseudocluster_E << " ) and " << second_highest_pseudocluster_index << " ( E = " << second_highest_pseudocluster_E << " ) " << std::endl;
1361  }
1362  // assign these clusters as owners
1363  tower_ownership[shared_ID] = std::pair<int, int>(highest_pseudocluster_index, second_highest_pseudocluster_index);
1364  }
1365 
1366  if (Verbosity() > 100)
1367  {
1368  for (int &original_tower : original_towers)
1369  {
1370  std::pair<int, int> the_pair = tower_ownership[original_tower];
1371  std::cout << " Debug Post-Split: tower_ownership[ " << original_tower << " ] = ( " << the_pair.first << ", " << the_pair.second << " ) ";
1372  std::cout << " , layer / ieta / iphi = " << get_ilayer_from_ID(original_tower) << " / " << get_ieta_from_ID(original_tower) << " / " << get_iphi_from_ID(original_tower);
1373  std::cout << std::endl;
1374  if (the_pair.first == -1)
1375  {
1376  std::vector<int> adjacent_tower_IDs = get_adjacent_towers_by_ID(original_tower);
1377 
1378  for (int this_adjacent_tower_ID : adjacent_tower_IDs)
1379  {
1380  if (get_status_from_ID(this_adjacent_tower_ID) != cl)
1381  {
1382  continue;
1383  }
1384  std::cout << " -> adjacent to add tower " << this_adjacent_tower_ID << " , which has status " << tower_ownership[this_adjacent_tower_ID].first << std::endl;
1385  }
1386  }
1387  }
1388  }
1389 
1390  // call helper function
1391  export_clusters(original_towers, tower_ownership, local_maxima_ID.size(), pseudocluster_sumE, pseudocluster_eta, pseudocluster_phi);
1392  }
1393 
1394  if (Verbosity() > 1)
1395  {
1396  std::cout << "RawClusterBuilderTopo::process_event after splitting (if any) final clusters output to node are: " << std::endl;
1397  RawClusterContainer::ConstRange begin_end = _clusters->getClusters();
1398  int ncl = 0;
1399  for (RawClusterContainer::ConstIterator hiter = begin_end.first; hiter != begin_end.second; ++hiter)
1400  {
1401  std::cout << "-> #" << ncl++ << " ";
1402  hiter->second->identify();
1403  std::cout << std::endl;
1404  }
1405  }
1406 
1407  return Fun4AllReturnCodes::EVENT_OK;
1408 }
1409 
1410 int RawClusterBuilderTopo::End(PHCompositeNode * /*topNode*/)
1411 {
1412  return Fun4AllReturnCodes::EVENT_OK;
1413 }
1414 
1415 void RawClusterBuilderTopo::CreateNodes(PHCompositeNode *topNode)
1416 {
1417  PHNodeIterator iter(topNode);
1418 
1419  PHCompositeNode *dstNode = static_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
1420  if (!dstNode)
1421  {
1422  std::cerr << PHWHERE << "DST Node missing, doing nothing." << std::endl;
1423  throw std::runtime_error("Failed to find DST node in RawClusterBuilderTopo::CreateNodes");
1424  }
1425 
1426  PHNodeIterator dstiter(dstNode);
1427  PHCompositeNode *DetNode = dynamic_cast<PHCompositeNode *>(dstiter.findFirst("PHCompositeNode", "TOPOCLUSTER"));
1428  if (!DetNode)
1429  {
1430  DetNode = new PHCompositeNode("TOPOCLUSTER");
1431  dstNode->addNode(DetNode);
1432  }
1433 
1434  _clusters = new RawClusterContainer();
1435 
1436  PHIODataNode<PHObject> *clusterNode = new PHIODataNode<PHObject>(_clusters, ClusterNodeName, "PHObject");
1437  DetNode->addNode(clusterNode);
1438 }