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