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>
36 return (a.second > b.second);
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 };
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 };
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 };
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 ) );
73 std::vector<int> adjacent_towers;
76 if ( this_layer == 0 || this_layer == 1 ) {
78 for (
int delta_layer = 0; delta_layer <= 1; delta_layer++) {
79 for (
int delta_eta = -1; delta_eta <= 1; delta_eta++) {
82 if ( delta_layer == 0 && delta_eta == 0 &&
delta_phi == 0 )
continue;
84 int test_eta = this_eta + delta_eta;
85 if ( test_eta < 0 || test_eta >=
_HCAL_NETA ) {
continue; }
87 int test_layer = ( this_layer + delta_layer ) % 2;
92 if (
Verbosity() > 20) std::cout <<
"RawClusterBuilderTopo::get_adjacent_towers_by_ID : corner growth not allowed " << std::endl;
97 adjacent_towers.push_back(
get_ID( test_layer, test_eta, test_phi ) );
111 for (
int new_eta = EMCal_eta_start; new_eta <= EMCal_eta_end; new_eta++) {
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 );
123 if ( this_layer == 2 ) {
126 for (
int delta_eta = -1; delta_eta <= 1; delta_eta++) {
129 if ( delta_eta == 0 &&
delta_phi == 0 )
continue;
131 int test_eta = this_eta + delta_eta;
132 if ( test_eta < 0 || test_eta >=
_EMCAL_NETA ) {
continue; }
137 adjacent_towers.push_back(
get_ID( this_layer, test_eta, test_phi ) );
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 );
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;
158 return adjacent_towers;
165 std::cout <<
"RawClusterBuilderTopo::export_single_cluster called " << std::endl;
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);
171 export_clusters( original_towers, tower_ownership, 1, std::vector<float>(), std::vector<float>(), std::vector<float>() );
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 ) {
179 if ( n_clusters != 1 )
181 std::cout <<
"RawClusterBuilderTopo::export_clusters called on an initial cluster with " << n_clusters <<
" final clusters " << std::endl;
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;
190 for (
unsigned int pc = 0;
pc < n_clusters;
pc++) {
192 clusters_E.push_back( 0 );
193 clusters_x.push_back( 0 );
194 clusters_y.push_back( 0 );
195 clusters_z.push_back( 0 );
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 ];
203 std::cout <<
"RawClusterBuilderTopo::export_clusters -> assigning tower " << original_towers.at(
t ) <<
" with ownership ( " << the_pair.first <<
", " << the_pair.second <<
" ) " << std::endl;
211 if ( this_layer == 2 ) {
220 if ( the_pair.second == -1 ) {
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();
229 std::cout <<
" -> tower ID " << this_ID <<
" fully assigned to pseudocluster " << the_pair.first << std::endl;
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 ] );
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;
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;
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);
260 for (
unsigned int cl = 0; cl < n_clusters; cl++) {
262 clusters[ cl ]->set_energy( clusters_E[ cl ] );
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 ];
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 );
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;
295 for (
int i=0; i<3; ++i)
327 catch (std::exception &
e)
329 std::cout <<
PHWHERE <<
": " << e.what() << std::endl;
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;
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");
353 std::cout <<
" RawClusterBuilderTopo::process_event : container TOWER_CALIB_CEMC does not exist, aborting " << std::endl;
357 std::cout <<
" RawClusterBuilderTopo::process_event : container TOWER_CALIB_HCALIN does not exist, aborting " << std::endl;
361 std::cout <<
" RawClusterBuilderTopo::process_event : container TOWER_CALIB_HCALOUT does not exist, aborting " << std::endl;
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");
369 if ( ! _geom_containers[0] ) {
370 std::cout <<
" RawClusterBuilderTopo::process_event : container TOWERGEOM_HCALIN does not exist, aborting " << std::endl;
373 if ( ! _geom_containers[1] ) {
374 std::cout <<
" RawClusterBuilderTopo::process_event : container TOWERGEOM_HCALOUT does not exist, aborting " << std::endl;
377 if ( !_geom_containers[2] ) {
378 std::cout <<
" RawClusterBuilderTopo::process_event : container TOWERGEOM_CEMC does not exist, aborting " << std::endl;
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;
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;
409 _HCAL_NETA = _geom_containers[1]->get_etabins();
410 _HCAL_NPHI = _geom_containers[1]->get_phibins();
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++) {
440 std::vector< std::pair<int, float> > list_of_seeds;
451 int ieta = _geom_containers[2]->get_etabin( tower_geom->
get_eta() );
452 int iphi = _geom_containers[2]->get_phibin( tower_geom->
get_phi() );
460 int ID =
get_ID( 2, ieta, iphi );
461 list_of_seeds.push_back( std::pair<int, float>( ID, this_E ) );
463 std::cout <<
"RawClusterBuilderTopo::process_event: adding EMCal tower at ieta / iphi = " << ieta <<
" / " << iphi <<
" with E = " << this_E << std::endl;
480 int ieta = _geom_containers[0]->get_etabin( tower_geom->
get_eta() );
481 int iphi = _geom_containers[0]->get_phibin( tower_geom->
get_phi() );
489 int ID =
get_ID( 0, ieta, iphi );
490 list_of_seeds.push_back( std::pair<int, float>( ID, this_E ) );
492 std::cout <<
"RawClusterBuilderTopo::process_event: adding IHCal tower at ieta / iphi = " << ieta <<
" / " << iphi <<
" with E = " << this_E << std::endl;
506 int ieta = _geom_containers[1]->get_etabin( tower_geom->
get_eta() );
507 int iphi = _geom_containers[1]->get_phibin( tower_geom->
get_phi() );
515 int ID =
get_ID( 1, ieta, iphi );
516 list_of_seeds.push_back( std::pair<int, float>( ID, this_E ) );
518 std::cout <<
"RawClusterBuilderTopo::process_event: adding OHCal tower at ieta / iphi = " << ieta <<
" / " << iphi <<
" with E = " << this_E << std::endl;
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;
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;
541 std::cout <<
"RawClusterBuilderTopo::process_event: initialized with " << list_of_seeds.size() <<
" seeds with E > 4*sigma " << std::endl;
544 int cluster_index = 0;
546 std::vector< std::vector<int> > all_cluster_towers;
548 while ( list_of_seeds.size() > 0 ) {
550 int seed_ID = list_of_seeds.at( 0 ).first;
551 list_of_seeds.erase( list_of_seeds.begin() );
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;
559 if ( seed_status > -1 ) {
561 std::cout <<
" --> already owned by cluster # " << seed_status << std::endl;
568 std::vector<int> cluster_tower_ID;
569 cluster_tower_ID.push_back( seed_ID );
571 std::vector<int> grow_tower_ID;
572 grow_tower_ID.push_back( seed_ID );
577 std::cout <<
" RawClusterBuilderTopo::process_event: Entering Growth stage for cluster " << cluster_index << std::endl;
579 while ( grow_tower_ID.size() > 0 ) {
581 int grow_ID = grow_tower_ID.at( 0 );
582 grow_tower_ID.erase( grow_tower_ID.begin() );
585 std::cout <<
" --> cluster " << cluster_index <<
", growth stage, examining neighbors of ID " << grow_ID <<
", " << grow_tower_ID.size() <<
" grow towers left" << std::endl;
589 for (
unsigned int adj = 0; adj < adjacent_tower_IDs.size(); adj++) {
591 int this_adjacent_tower_ID = adjacent_tower_IDs.at( adj );
593 if (
Verbosity() > 10) std::cout <<
" --> --> --> checking possible adjacent tower with ID " << this_adjacent_tower_ID <<
" : ";
599 if (
Verbosity() > 10) std::cout <<
"does not exist " << std::endl;
605 if (
Verbosity() > 10) std::cout <<
"already owned by this cluster index " << cluster_index << std::endl;
611 if (
Verbosity() > 10) std::cout <<
"E = " <<
get_E_from_ID( this_adjacent_tower_ID ) <<
" under 2*sigma threshold " << std::endl;
617 if (
Verbosity() > 10) std::cout <<
"ERROR! in growth stage, encountered >2sigma tower which is already owned?!" << std::endl;
622 grow_tower_ID.push_back( this_adjacent_tower_ID );
623 cluster_tower_ID.push_back( this_adjacent_tower_ID );
625 if (
Verbosity() > 10) std::cout <<
"add this tower ( ID " << this_adjacent_tower_ID <<
" ) to grow list " << std::endl;
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;
635 std::cout <<
" RawClusterBuilderTopo::process_event: Entering Perimeter stage for cluster " << cluster_index << std::endl;
638 int n_core_towers = cluster_tower_ID.size();
640 for (
int ic = 0; ic < n_core_towers; ic++) {
642 int core_ID = cluster_tower_ID.at( ic );
645 std::cout <<
" --> cluster " << cluster_index <<
", perimeter stage, examining neighbors of ID " << core_ID <<
", core cluster # " << ic <<
" of " << n_core_towers <<
" total " << std::endl;
649 for (
unsigned int adj = 0; adj < adjacent_tower_IDs.size(); adj++) {
651 int this_adjacent_tower_ID = adjacent_tower_IDs.at( adj );
653 if (
Verbosity() > 10) std::cout <<
" --> --> --> checking possible adjacent tower with ID " << this_adjacent_tower_ID <<
" : ";
659 if (
Verbosity() > 10) std::cout <<
"does not exist " << std::endl;
665 if (
Verbosity() > 10) std::cout <<
"already owned by other cluster index " <<
get_status_from_ID( this_adjacent_tower_ID ) << std::endl;
671 if (
Verbosity() > 10) std::cout <<
"E = " <<
get_E_from_ID( this_adjacent_tower_ID ) <<
" under 0*sigma threshold " << std::endl;
676 cluster_tower_ID.push_back( this_adjacent_tower_ID );
678 if (
Verbosity() > 10) std::cout <<
"add this tower ( ID " << this_adjacent_tower_ID <<
" ) to cluster " << std::endl;
682 if (
Verbosity() > 5) std::cout <<
" --> after examining perimeter neighbors, # of towers in cluster is now = " << cluster_tower_ID.size() << std::endl;
686 all_cluster_towers.push_back( cluster_tower_ID );
693 if (
Verbosity() > 0) std::cout <<
"RawClusterBuilderTopo::process_event: " << cluster_index <<
" topo-clusters initially reconstructed, entering splitting step" << std::endl;
696 int original_cluster_index = cluster_index;
697 for (
int cl = 0; cl < original_cluster_index; cl++) {
699 std::vector<int> original_towers = all_cluster_towers.at( cl );
703 if (
Verbosity() > 2 ) std::cout <<
"RawClusterBuilderTopo::process_event: splitting step disabled, cluster " << cluster_index <<
" is final" << std::endl;
708 std::vector< std::pair<int, float> > local_maxima_ID;
711 for (
unsigned int t = 0;
t < original_towers.size();
t++) {
712 int tower_ID = original_towers.at(
t );
714 if (
Verbosity() > 10 ) std::cout <<
" -> examining tower ID " << tower_ID <<
" for possible local maximum " << std::endl;
724 int neighbors_in_cluster = 0;
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 );
733 neighbors_in_cluster++;
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;
742 if (has_higher_neighbor)
continue;
745 if ( neighbors_in_cluster < 4 ) {
746 if (
Verbosity() > 10 ) std::cout <<
" -> -> too few neighbors N = " << neighbors_in_cluster << std::endl;
750 local_maxima_ID.push_back( std::pair<int,float>( tower_ID ,
get_E_from_ID( tower_ID ) ) );
756 for (
unsigned int n = 0;
n < local_maxima_ID.size();
n++) {
759 std::pair<int,float> this_LM = local_maxima_ID.at(
n );
763 if ( this_phi > 3.14159 ) this_phi -= 2 * 3.14159;
766 bool has_EM_overlap =
false;
769 for (
unsigned int n2 = 0; n2 < local_maxima_ID.size(); n2++) {
771 if (
n == n2 )
continue;
774 std::pair<int,float> this_LM2 = local_maxima_ID.at( n2 );
778 if ( this_phi2 > 3.14159 ) this_phi -= 2 * 3.14159;
782 float dR =
calculate_dR( this_eta, this_eta2, this_phi, this_phi2 );
786 has_EM_overlap =
true;
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;
795 if ( has_EM_overlap ) {
797 local_maxima_ID.erase( local_maxima_ID.begin() +
n );
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 ) <<
", ";
810 if ( this_phi > 3.14159 ) this_phi -= 2 * 3.14159;
817 if ( local_maxima_ID.size() <= 1 ) {
819 if (
Verbosity() > 2) std::cout <<
"RawClusterBuilderTopo::process_event cluster " << cl <<
" has only " << local_maxima_ID.size() <<
" local maxima, not splitting " << std::endl;
829 std::cout <<
"RawClusterBuilderTopo::process_event splitting cluster " << cl <<
" into " << local_maxima_ID.size() <<
" according to local maxima!" << std::endl;
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);
839 std::vector<int> seed_list;
840 std::vector<int> neighbor_list;
841 std::vector<int> shared_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 );
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 <<
" ) ";
858 std::cout << std::endl;
862 bool first_pass =
true;
867 std::cout <<
" -> starting split loop with " << seed_list.size() <<
" seed, " << neighbor_list.size() <<
" neighbor, and " << shared_list.size() <<
" shared towers " << std::endl;
870 std::vector<int> new_ownerships;
872 for (
unsigned int n = 0;
n < neighbor_list.size();
n++) {
874 int neighbor_ID = neighbor_list.at(
n );
877 std::cout <<
" -> -> looking at neighbor " <<
n <<
" (tower ID " << neighbor_ID <<
" ) of " << neighbor_list.size() <<
" total" << std::endl;
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 );
885 std::map<int, bool> pseudocluster_adjacency;
886 for (
unsigned int s = 0;
s < local_maxima_ID.size();
s++)
887 pseudocluster_adjacency[
s ] =
false;
891 for (
unsigned int adj = 0; adj < adjacent_tower_IDs.size(); adj++) {
892 int this_adjacent_tower_ID = adjacent_tower_IDs.at( adj );
894 if ( tower_ownership[ this_adjacent_tower_ID ].first > -1 ) {
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;
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++;
907 std::cout <<
" -> -> adjacent to pseudocluster " << s << std::endl;
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 );
915 else if ( n_pseudocluster_adjacent == 1 ) {
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 );
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 );
929 std::cout <<
" -> now updating status of all " << neighbor_list.size() <<
" original neighbors " << std::endl;
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 );
937 std::cout <<
" -> -> neighbor ID " << neighbor_ID <<
" has new status " << new_ownerships.at(
n ) << std::endl;
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 );
943 std::cout <<
" -> -> neighbor ID " << neighbor_ID <<
" has new status " << -3 << std::endl;
948 std::cout <<
" producing a new neighbor list ... " << std::endl;
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 ) {
956 for (
unsigned int adj = 0; adj < adjacent_tower_IDs.size(); adj++) {
957 int this_adjacent_tower_ID = adjacent_tower_IDs.at( adj );
959 if ( tower_ownership[ this_adjacent_tower_ID ].first == -1 ) {
960 new_neighbor_list.push_back( this_adjacent_tower_ID );
962 std::cout <<
" -> queueing up to add tower " << this_adjacent_tower_ID <<
" , neighbor of tower " << neighbor_ID <<
" to new neighbor list" << std::endl;
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;
977 neighbor_list.clear();
980 for (; new_neighbor_list.size() > 0; ) {
981 neighbor_list.push_back( new_neighbor_list.front( ) );
982 new_neighbor_list.pop_front();
987 }
while ( neighbor_list.size() > 0 );
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 <<
" ) ";
994 std::cout << std::endl;
995 if ( the_pair.first == -1 ) {
998 for (
unsigned int adj = 0; adj < adjacent_tower_IDs.size(); adj++) {
999 int this_adjacent_tower_ID = adjacent_tower_IDs.at( adj );
1001 std::cout <<
" -> adjacent to add tower " << this_adjacent_tower_ID <<
" , which has status " << tower_ownership[ this_adjacent_tower_ID ].first << std::endl;
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;
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 );
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 );
1029 pseudocluster_sumeta[ the_pair.first ] += this_eta;
1030 pseudocluster_sumphi[ the_pair.first ] += this_phi;
1031 pseudocluster_ntower[ the_pair.first ] += 1;
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 ) );
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;
1045 std::cout <<
"RawClusterBuilderTopo::process_event now splitting up shared clusters (including unassigned clusters), initial shared list has size " << shared_list.size() << std::endl;
1048 while ( shared_list.size() > 0 ) {
1051 int shared_ID = shared_list.at( 0 );
1052 shared_list.erase( shared_list.begin() );
1055 std::cout <<
" -> looking at shared tower " << shared_ID <<
", after this one there are " << shared_list.size() <<
" shared towers left " << std::endl;
1058 std::vector<bool> pseudocluster_adjacency;
1059 pseudocluster_adjacency.resize( local_maxima_ID.size(),
false );
1063 for (
unsigned int adj = 0; adj < adjacent_tower_IDs.size(); adj++) {
1064 int this_adjacent_tower_ID = adjacent_tower_IDs.at( adj );
1066 if ( tower_ownership[ this_adjacent_tower_ID ].first > -1 ) {
1067 pseudocluster_adjacency[ tower_ownership[ this_adjacent_tower_ID ].first ] =
true;
1069 if ( tower_ownership[ this_adjacent_tower_ID ].
second > -1 ) {
1070 pseudocluster_adjacency[ tower_ownership[ this_adjacent_tower_ID ].second ] =
true;
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);
1077 std::cout <<
" -> while looking at neighbors, have added un-examined tower " << this_adjacent_tower_ID <<
" to shared list " << std::endl;
1082 int highest_pseudocluster_index = -1;
1083 int second_highest_pseudocluster_index = -1;
1085 float highest_pseudocluster_E = -1;
1086 float second_highest_pseudocluster_E = -2;
1088 for (
unsigned int n = 0;
n < pseudocluster_adjacency.size();
n++) {
1090 if ( ! pseudocluster_adjacency[
n ] )
continue;
1092 if ( pseudocluster_sumE[ n ] > highest_pseudocluster_E ) {
1093 second_highest_pseudocluster_E = highest_pseudocluster_E;
1094 second_highest_pseudocluster_index = highest_pseudocluster_index;
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;
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;
1109 tower_ownership[ shared_ID ] = std::pair<int, int>( highest_pseudocluster_index, second_highest_pseudocluster_index );
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 <<
" ) ";
1118 std::cout << std::endl;
1119 if ( the_pair.first == -1 ) {
1122 for (
unsigned int adj = 0; adj < adjacent_tower_IDs.size(); adj++) {
1123 int this_adjacent_tower_ID = adjacent_tower_IDs.at( adj );
1125 std::cout <<
" -> adjacent to add tower " << this_adjacent_tower_ID <<
" , which has status " << tower_ownership[ this_adjacent_tower_ID ].first << std::endl;
1133 export_clusters( original_towers , tower_ownership , local_maxima_ID.size() , pseudocluster_sumE, pseudocluster_eta, pseudocluster_phi );
1138 std::cout <<
"RawClusterBuilderTopo::process_event after splitting (if any) final clusters output to node are: " << std::endl;
1143 std::cout <<
"-> #" << ncl++ <<
" " ;
1144 hiter->second->identify();
1145 std::cout << std::endl;
1164 std::cerr <<
PHWHERE <<
"DST Node missing, doing nothing." << std::endl;
1165 throw std::runtime_error(
"Failed to find DST node in RawClusterBuilderTopo::CreateNodes");
1179 DetNode->
addNode(clusterNode);