Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-12-17 09:19:59 

0001 #include "RawClusterBuilderTopo.h"
0002 
0003 #include <calobase/RawCluster.h>
0004 #include <calobase/RawClusterContainer.h>
0005 #include <calobase/RawClusterv1.h>
0006 #include <calobase/RawTowerDefs.h>  // for encode_towerid, Calorime...
0007 #include <calobase/RawTowerGeom.h>
0008 #include <calobase/RawTowerGeomContainer.h>
0009 #include <calobase/TowerInfo.h>           // for TowerInfo
0010 #include <calobase/TowerInfoContainer.h>  // for TowerInfoContainer
0011 
0012 #include <fun4all/Fun4AllReturnCodes.h>
0013 #include <fun4all/SubsysReco.h>
0014 
0015 #include <phool/PHCompositeNode.h>
0016 #include <phool/PHIODataNode.h>
0017 #include <phool/PHNode.h>
0018 #include <phool/PHNodeIterator.h>
0019 #include <phool/PHObject.h>
0020 #include <phool/getClass.h>
0021 #include <phool/phool.h>
0022 
0023 #include <algorithm>
0024 #include <cmath>
0025 #include <cstdlib>  // for abs
0026 #include <exception>
0027 #include <iostream>
0028 #include <iterator>  // for begin, end
0029 #include <list>
0030 #include <memory>  // for allocator_traits<>::valu...
0031 #include <stdexcept>
0032 #include <utility>
0033 #include <vector>
0034 
0035 namespace
0036 {
0037 bool sort_by_pair_second(const std::pair<int, float> &a, const std::pair<int, float> &b)
0038 {
0039   return (a.second > b.second);
0040 }
0041 }
0042 
0043 int RawClusterBuilderTopo::RawClusterBuilderTopo_constants_EMCal_eta_start_given_IHCal[24] =
0044     {2, 6, 10, 14, 18, 22, 26, 30, 33, 37, 41, 44,
0045      48, 52, 55, 59, 63, 66, 70, 74, 78, 82, 86, 90};
0046 
0047 int RawClusterBuilderTopo::RawClusterBuilderTopo_constants_EMCal_eta_end_given_IHCal[24] =
0048     {5, 9, 13, 17, 21, 25, 29, 32, 36, 40, 43, 47,
0049      51, 54, 58, 62, 65, 69, 73, 77, 81, 85, 89, 93};
0050 
0051 int RawClusterBuilderTopo::RawClusterBuilderTopo_constants_IHCal_eta_given_EMCal[96] = {
0052     -1, -1, 0, 0, 0, 0, 1, 1, 1, 1, 2, 2,
0053     2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 5, 5,
0054     5, 5, 6, 6, 6, 6, 7, 7, 7, 8, 8, 8,
0055     8, 9, 9, 9, 9, 10, 10, 10, 11, 11, 11, 11,
0056     12, 12, 12, 12, 13, 13, 13, 14, 14, 14, 14, 15,
0057     15, 15, 15, 16, 16, 16, 17, 17, 17, 17, 18, 18,
0058     18, 18, 19, 19, 19, 19, 20, 20, 20, 20, 21, 21,
0059     21, 21, 22, 22, 22, 22, 23, 23, 23, 23, -1, -1};
0060 
0061 float RawClusterBuilderTopo::calculate_dR(float eta1, float eta2, float phi1, float phi2)
0062 {
0063   float deta = eta1 - eta2;
0064   float dphi = phi1 - phi2;
0065   while (dphi > M_PI)
0066   {
0067     dphi -= 2 * M_PI;
0068   }
0069   while (dphi < -M_PI)
0070   {
0071     dphi += 2 * M_PI;
0072   }
0073   return sqrt(pow(deta, 2) + pow(dphi, 2));
0074 }
0075 
0076 std::vector<int> RawClusterBuilderTopo::get_adjacent_towers_by_ID(int ID)
0077 {
0078   int this_layer = get_ilayer_from_ID(ID);
0079   int this_eta = get_ieta_from_ID(ID);
0080   int this_phi = get_iphi_from_ID(ID);
0081 
0082   std::vector<int> adjacent_towers;
0083 
0084   // for both IHCal and OHCal, add adjacent layers in the HCal
0085   if (this_layer == 0 || this_layer == 1)
0086   {
0087     for (int delta_layer = 0; delta_layer <= 1; delta_layer++)
0088     {
0089       for (int delta_eta = -1; delta_eta <= 1; delta_eta++)
0090       {
0091         for (int delta_phi = -1; delta_phi <= 1; delta_phi++)
0092         {
0093           if (delta_layer == 0 && delta_eta == 0 && delta_phi == 0)
0094           {
0095             continue;  // this is the same tower
0096           }
0097 
0098           int test_eta = this_eta + delta_eta;
0099           if (test_eta < 0 || test_eta >= _HCAL_NETA)
0100           {
0101             continue;
0102           }  // ignore if at the (eta) edge of calorimeter
0103 
0104           int test_layer = (this_layer + delta_layer) % 2;                  // wrap around in layer
0105           int test_phi = (this_phi + delta_phi + _HCAL_NPHI) % _HCAL_NPHI;  // wrap around in phi (add 64 to avoid -1)
0106 
0107           // disallow "corner" adjacency (diagonal in eta/phi plane and in different layer) if this option not enabled
0108           if (!_allow_corner_neighbor && delta_layer == 1 && abs(delta_phi) == 1 && abs(delta_eta) == 1)
0109           {
0110             if (Verbosity() > 20)
0111             {
0112               std::cout << "RawClusterBuilderTopo::get_adjacent_towers_by_ID : corner growth not allowed " << std::endl;
0113             }
0114             continue;
0115           }
0116 
0117           // add to list of adjacent towers
0118           adjacent_towers.push_back(get_ID(test_layer, test_eta, test_phi));
0119         }
0120       }
0121     }
0122   }
0123 
0124   // for IHCal only, also add 4x4 group of EMCal towers
0125   if (this_layer == 0 && _enable_EMCal)
0126   {
0127     int EMCal_phi_start = get_first_matching_EMCal_phi_from_IHCal(this_phi);
0128     int EMCal_eta_start = RawClusterBuilderTopo_constants_EMCal_eta_start_given_IHCal[this_eta];
0129     int EMCal_eta_end = RawClusterBuilderTopo_constants_EMCal_eta_end_given_IHCal[this_eta];
0130 
0131     for (int new_eta = EMCal_eta_start; new_eta <= EMCal_eta_end; new_eta++)
0132     {
0133       for (int delta_phi = 0; delta_phi < 4; delta_phi++)
0134       {
0135         int new_phi = (EMCal_phi_start + delta_phi + _EMCAL_NPHI) % _EMCAL_NPHI;
0136 
0137         int EMCal_tower = get_ID(2, new_eta, new_phi);
0138         if (Verbosity() > 20)
0139         {
0140           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;
0141         }
0142         adjacent_towers.push_back(EMCal_tower);
0143       }
0144     }
0145   }
0146 
0147   // for EMCal, add adjacent EMCal towers and (usually) one IHCal tower
0148   if (this_layer == 2)
0149   {
0150     // EMCal towers first
0151     for (int delta_eta = -1; delta_eta <= 1; delta_eta++)
0152     {
0153       for (int delta_phi = -1; delta_phi <= 1; delta_phi++)
0154       {
0155         if (delta_eta == 0 && delta_phi == 0)
0156         {
0157           continue;  // this is the same tower
0158         }
0159 
0160         int test_eta = this_eta + delta_eta;
0161         if (test_eta < 0 || test_eta >= _EMCAL_NETA)
0162         {
0163           continue;
0164         }  // ignore if at the (eta) edge of calorimeter
0165 
0166         int test_phi = (this_phi + delta_phi + _EMCAL_NPHI) % _EMCAL_NPHI;  // wrap around in phi (add 256 to avoid -1)
0167 
0168         // add to list of adjacent towers
0169         adjacent_towers.push_back(get_ID(this_layer, test_eta, test_phi));
0170       }
0171     }
0172 
0173     // now add IHCal towers
0174     if (_enable_HCal)
0175     {
0176       int HCal_eta = RawClusterBuilderTopo_constants_IHCal_eta_given_EMCal[this_eta];
0177       int HCal_phi = get_matching_HCal_phi_from_EMCal(this_phi);
0178 
0179       if (HCal_eta >= 0)
0180       {
0181         int IHCal_tower = get_ID(0, HCal_eta, HCal_phi);
0182         if (Verbosity() > 20)
0183         {
0184           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;
0185         }
0186         adjacent_towers.push_back(IHCal_tower);
0187       }
0188       else
0189       {
0190         if (Verbosity() > 20)
0191         {
0192           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;
0193         }
0194       }
0195     }
0196   }
0197 
0198   return adjacent_towers;
0199 }
0200 
0201 void RawClusterBuilderTopo::export_single_cluster(const std::vector<int> &original_towers)
0202 {
0203   if (Verbosity() > 2)
0204   {
0205     std::cout << "RawClusterBuilderTopo::export_single_cluster called " << std::endl;
0206   }
0207 
0208   std::map<int, std::pair<int, int> > tower_ownership;
0209   for (const int &original_tower : original_towers)
0210   {
0211     tower_ownership[original_tower] = std::pair<int, int>(0, -1);  // all towers owned by cluster 0
0212   }
0213   export_clusters(original_towers, tower_ownership, 1, std::vector<float>(), std::vector<float>(), std::vector<float>());
0214 
0215   return;
0216 }
0217 
0218 void RawClusterBuilderTopo::export_clusters(const std::vector<int> &original_towers, std::map<int, std::pair<int, int> > tower_ownership, unsigned int n_clusters, const std::vector<float> &pseudocluster_sumE, const std::vector<float> &pseudocluster_eta, const std::vector<float> &pseudocluster_phi)
0219 {
0220   if (n_clusters != 1)  // if we didn't just pass down from export_single_cluster
0221   {
0222     if (Verbosity() > 2)
0223     {
0224       std::cout << "RawClusterBuilderTopo::export_clusters called on an initial cluster with " << n_clusters << " final clusters " << std::endl;
0225     }
0226   }
0227   // build a RawCluster for output
0228   std::vector<RawCluster *> clusters;
0229   std::vector<float> clusters_E;
0230   std::vector<float> clusters_absE;
0231   std::vector<float> clusters_x;
0232   std::vector<float> clusters_y;
0233   std::vector<float> clusters_z;
0234 
0235   for (unsigned int pc = 0; pc < n_clusters; pc++)
0236   {
0237     clusters.push_back(new RawClusterv1());
0238     clusters_E.push_back(0);
0239     clusters_absE.push_back(0);
0240     clusters_x.push_back(0);
0241     clusters_y.push_back(0);
0242     clusters_z.push_back(0);
0243   }
0244 
0245   for (int original_tower : original_towers)
0246   {
0247     int this_ID = original_tower;
0248     std::pair<int, int> the_pair = tower_ownership[this_ID];
0249 
0250     if (Verbosity() > 5)
0251     {
0252       std::cout << "RawClusterBuilderTopo::export_clusters -> assigning tower " << original_tower << " with ownership ( " << the_pair.first << ", " << the_pair.second << " ) " << std::endl;
0253     }
0254     int this_ieta = get_ieta_from_ID(this_ID);
0255     int this_iphi = get_iphi_from_ID(this_ID);
0256     int this_layer = get_ilayer_from_ID(this_ID);
0257     float this_E = get_E_from_ID(this_ID);
0258 
0259     int this_key = 0;
0260     if (this_layer == 2)
0261     {
0262       this_key = _EMTOWERMAP_KEY_ETA_PHI[this_ieta][this_iphi];
0263     }
0264     else
0265     {
0266       this_key = _TOWERMAP_KEY_LAYER_ETA_PHI[this_layer][this_ieta][this_iphi];
0267     }
0268 
0269     RawTowerGeom *tower_geom = _geom_containers[this_layer]->get_tower_geometry(this_key);
0270 
0271     if (the_pair.second == -1)
0272     {
0273       // assigned only to one cluster, easy
0274       clusters[the_pair.first]->addTower(this_key, this_E);
0275       clusters_E[the_pair.first] = clusters_E[the_pair.first] + this_E;
0276       clusters_absE[the_pair.first] = clusters_absE[the_pair.first] + std::fabs(this_E);
0277       // calculate position mean using absolute energy as weights
0278       clusters_x[the_pair.first] = clusters_x[the_pair.first] + std::fabs(this_E) * tower_geom->get_center_x();
0279       clusters_y[the_pair.first] = clusters_y[the_pair.first] + std::fabs(this_E) * tower_geom->get_center_y();
0280       clusters_z[the_pair.first] = clusters_z[the_pair.first] + std::fabs(this_E) * tower_geom->get_center_z();
0281 
0282       if (Verbosity() > 5)
0283       {
0284         std::cout << " -> tower ID " << this_ID << " fully assigned to pseudocluster " << the_pair.first << std::endl;
0285       }
0286     }
0287     else
0288     {
0289       // assigned to two clusters! get energy sharing fraction ...
0290       float dR1 = calculate_dR(tower_geom->get_eta(), pseudocluster_eta[the_pair.first], tower_geom->get_phi(), pseudocluster_phi[the_pair.first]) / _R_shower;
0291       float dR2 = calculate_dR(tower_geom->get_eta(), pseudocluster_eta[the_pair.second], tower_geom->get_phi(), pseudocluster_phi[the_pair.second]) / _R_shower;
0292       float r = std::exp(dR1 - dR2);
0293       float frac1 = fabs(pseudocluster_sumE[the_pair.first]) / (fabs(pseudocluster_sumE[the_pair.first]) + r * fabs(pseudocluster_sumE[the_pair.second]));
0294 
0295       if (Verbosity() > 5)
0296       {
0297         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;
0298       }
0299       clusters[the_pair.first]->addTower(this_key, this_E * frac1);
0300       clusters_E[the_pair.first] = clusters_E[the_pair.first] + this_E * frac1;
0301       clusters_absE[the_pair.first] = clusters_absE[the_pair.first] + std::fabs(this_E) * frac1;
0302       clusters_x[the_pair.first] = clusters_x[the_pair.first] + std::fabs(this_E) * tower_geom->get_center_x() * frac1;
0303       clusters_y[the_pair.first] = clusters_y[the_pair.first] + std::fabs(this_E) * tower_geom->get_center_y() * frac1;
0304       clusters_z[the_pair.first] = clusters_z[the_pair.first] + std::fabs(this_E) * tower_geom->get_center_z() * frac1;
0305 
0306       clusters[the_pair.second]->addTower(this_key, this_E * (1 - frac1));
0307       clusters_E[the_pair.second] = clusters_E[the_pair.second] + this_E * (1 - frac1);
0308       clusters_absE[the_pair.second] = clusters_absE[the_pair.second] + std::fabs(this_E) * (1 - frac1);
0309       clusters_x[the_pair.second] = clusters_x[the_pair.second] + std::fabs(this_E) * tower_geom->get_center_x() * (1 - frac1);
0310       clusters_y[the_pair.second] = clusters_y[the_pair.second] + std::fabs(this_E) * tower_geom->get_center_y() * (1 - frac1);
0311       clusters_z[the_pair.second] = clusters_z[the_pair.second] + std::fabs(this_E) * tower_geom->get_center_z() * (1 - frac1);
0312     }
0313   }
0314 
0315   // iterate through and add to official container
0316 
0317   for (unsigned int cl = 0; cl < n_clusters; cl++)
0318   {
0319     if (clusters_absE[cl] < _min_cluster_E)
0320     {
0321       if (Verbosity() > 2)
0322       {
0323         std::cout << "RawClusterBuilderTopo::export_clusters: skipping cluster with E = " << clusters_E[cl] << " and absE = " << clusters_absE[cl] << " due to low energy " << std::endl;
0324       }
0325       continue;
0326     }
0327     clusters[cl]->set_energy(clusters_E[cl]);
0328 
0329     float mean_x = clusters_x[cl] / clusters_absE[cl];
0330     float mean_y = clusters_y[cl] / clusters_absE[cl];
0331     float mean_z = clusters_z[cl] / clusters_absE[cl];
0332 
0333     clusters[cl]->set_r(std::sqrt((mean_y * mean_y) + (mean_x * mean_x)));
0334     clusters[cl]->set_phi(std::atan2(mean_y, mean_x));
0335     clusters[cl]->set_z(mean_z);
0336 
0337     _clusters->AddCluster(clusters[cl]);
0338 
0339     if (Verbosity() > 1)
0340     {
0341       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;
0342     }
0343   }
0344 
0345   return;
0346 }
0347 
0348 RawClusterBuilderTopo::RawClusterBuilderTopo(const std::string &name)
0349   : SubsysReco(name)
0350 {
0351   // geometry defined at run-time
0352 
0353   std::fill(std::begin(_geom_containers), std::end(_geom_containers), nullptr);
0354   _noise_LAYER[0] = 0.0025;
0355   _noise_LAYER[1] = 0.006;
0356   _noise_LAYER[2] = 0.03;  // EM
0357 
0358   _local_max_minE_LAYER[0] = 1;
0359   _local_max_minE_LAYER[1] = 1;
0360   _local_max_minE_LAYER[2] = 1;
0361 }
0362 
0363 int RawClusterBuilderTopo::InitRun(PHCompositeNode *topNode)
0364 {
0365   try
0366   {
0367     CreateNodes(topNode);
0368   }
0369   catch (std::exception &e)
0370   {
0371     std::cout << PHWHERE << ": " << e.what() << std::endl;
0372     throw;
0373   }
0374 
0375   if (Verbosity() > 0)
0376   {
0377     std::cout << "RawClusterBuilderTopo::InitRun: initialized with EMCal enable = " << _enable_EMCal << " and I+OHCal enable = " << _enable_HCal << std::endl;
0378     std::cout << "RawClusterBuilderTopo::InitRun: initialized with sigma_noise in EMCal / IHCal / OHCal = " << _noise_LAYER[2] << " / " << _noise_LAYER[0] << " / " << _noise_LAYER[1] << std::endl;
0379     std::cout << "RawClusterBuilderTopo::InitRun: initialized with noise multiples for seeding / growth / perimeter ( S / N / P ) = " << _sigma_seed << " / " << _sigma_grow << " / " << _sigma_peri << std::endl;
0380     std::cout << "RawClusterBuilderTopo::InitRun: initialized with allow_corner_neighbor = " << _allow_corner_neighbor << " (in HCal)" << std::endl;
0381     std::cout << "RawClusterBuilderTopo::InitRun: initialized with use_absE = " << _use_absE << std::endl;
0382     std::cout << "RawClusterBuilderTopo::InitRun: initialized with do_split = " << _do_split << " , R_shower = " << _R_shower << " (angular units) " << std::endl;
0383     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;
0384   }
0385 
0386   return Fun4AllReturnCodes::EVENT_OK;
0387 }
0388 
0389 int RawClusterBuilderTopo::process_event(PHCompositeNode *topNode)
0390 {
0391 
0392   std::string towerinfoNodenameEM = "TOWERINFO_CALIB_CEMC";
0393   std::string towerinfoNodenameIH = "TOWERINFO_CALIB_HCALIN";
0394   std::string towerinfoNodenameOH = "TOWERINFO_CALIB_HCALOUT";
0395   if (!_inputnodeprefix.empty())
0396   {
0397     towerinfoNodenameEM = _inputnodeprefix + "_CEMC";
0398     towerinfoNodenameIH = _inputnodeprefix + "_HCALIN";
0399     towerinfoNodenameOH = _inputnodeprefix + "_HCALOUT";
0400   }
0401 
0402   TowerInfoContainer *towerinfosEM = findNode::getClass<TowerInfoContainer>(topNode, towerinfoNodenameEM);
0403   TowerInfoContainer *towerinfosIH = findNode::getClass<TowerInfoContainer>(topNode, towerinfoNodenameIH);
0404   TowerInfoContainer *towerinfosOH = findNode::getClass<TowerInfoContainer>(topNode, towerinfoNodenameOH);
0405 
0406   if (!towerinfosEM)
0407   {
0408     std::cout << " RawClusterBuilderTopo::process_event : container TOWERINFO_CALIB_CEMC does not exist, aborting " << std::endl;
0409     return Fun4AllReturnCodes::ABORTEVENT;
0410   }
0411   if (!towerinfosIH)
0412   {
0413     std::cout << " RawClusterBuilderTopo::process_event : container TOWERINFO_CALIB_HCALIN does not exist, aborting " << std::endl;
0414     return Fun4AllReturnCodes::ABORTEVENT;
0415   }
0416   if (!towerinfosOH)
0417   {
0418     std::cout << " RawClusterBuilderTopo::process_event : container TOWERINFO_CALIB_HCALOUT does not exist, aborting " << std::endl;
0419     return Fun4AllReturnCodes::ABORTEVENT;
0420   }
0421 
0422   _geom_containers[0] = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALIN");
0423   _geom_containers[1] = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALOUT");
0424   _geom_containers[2] = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_CEMC");
0425 
0426   if (!_geom_containers[0])
0427   {
0428     std::cout << " RawClusterBuilderTopo::process_event : container TOWERGEOM_HCALIN does not exist, aborting " << std::endl;
0429     return Fun4AllReturnCodes::ABORTEVENT;
0430   }
0431   if (!_geom_containers[1])
0432   {
0433     std::cout << " RawClusterBuilderTopo::process_event : container TOWERGEOM_HCALOUT does not exist, aborting " << std::endl;
0434     return Fun4AllReturnCodes::ABORTEVENT;
0435   }
0436   if (!_geom_containers[2])
0437   {
0438     std::cout << " RawClusterBuilderTopo::process_event : container TOWERGEOM_CEMC does not exist, aborting " << std::endl;
0439     return Fun4AllReturnCodes::ABORTEVENT;
0440   }
0441 
0442   if (Verbosity() > 10)
0443   {
0444     std::cout << "RawClusterBuilderTopo::process_event: " << towerinfosEM->size() << " TOWERINFO_CALIB_CEMC towers" << std::endl;
0445     std::cout << "RawClusterBuilderTopo::process_event: " << towerinfosIH->size() << " TOWERINFO_CALIB_HCALIN towers" << std::endl;
0446     std::cout << "RawClusterBuilderTopo::process_event: " << towerinfosOH->size() << " TOWERINFO_CALIB_HCALOUT towers" << std::endl;
0447 
0448     std::cout << "RawClusterBuilderTopo::process_event: pointer to TOWERGEOM_CEMC: " << _geom_containers[2] << std::endl;
0449     std::cout << "RawClusterBuilderTopo::process_event: pointer to TOWERGEOM_HCALIN: " << _geom_containers[0] << std::endl;
0450     std::cout << "RawClusterBuilderTopo::process_event: pointer to TOWERGEOM_HCALOUT: " << _geom_containers[1] << std::endl;
0451   }
0452 
0453   if (_EMCAL_NETA < 0)
0454   {
0455     // define geometry only once if it has not been yet
0456     _EMCAL_NETA = _geom_containers[2]->get_etabins();
0457     _EMCAL_NPHI = _geom_containers[2]->get_phibins();
0458 
0459     _EMTOWERMAP_STATUS_ETA_PHI.resize(_EMCAL_NETA, std::vector<int>(_EMCAL_NPHI, -2));
0460     _EMTOWERMAP_KEY_ETA_PHI.resize(_EMCAL_NETA, std::vector<int>(_EMCAL_NPHI, 0));
0461     _EMTOWERMAP_E_ETA_PHI.resize(_EMCAL_NETA, std::vector<float>(_EMCAL_NPHI, 0));
0462   }
0463 
0464   if (_HCAL_NETA < 0)
0465   {
0466     // define geometry only once if it has not been yet
0467     _HCAL_NETA = _geom_containers[1]->get_etabins();
0468     _HCAL_NPHI = _geom_containers[1]->get_phibins();
0469 
0470     _TOWERMAP_STATUS_LAYER_ETA_PHI.resize(2, std::vector<std::vector<int> >(_HCAL_NETA, std::vector<int>(_HCAL_NPHI, -2)));
0471     _TOWERMAP_KEY_LAYER_ETA_PHI.resize(2, std::vector<std::vector<int> >(_HCAL_NETA, std::vector<int>(_HCAL_NPHI, 0)));
0472     _TOWERMAP_E_LAYER_ETA_PHI.resize(2, std::vector<std::vector<float> >(_HCAL_NETA, std::vector<float>(_HCAL_NPHI, 0)));
0473   }
0474 
0475   // reset maps
0476   // but note -- do not reset keys!
0477   for (int ieta = 0; ieta < _EMCAL_NETA; ieta++)
0478   {
0479     for (int iphi = 0; iphi < _EMCAL_NPHI; iphi++)
0480     {
0481       _EMTOWERMAP_STATUS_ETA_PHI[ieta][iphi] = -2;  // set tower does not exist
0482       _EMTOWERMAP_E_ETA_PHI[ieta][iphi] = 0;        // set zero energy
0483     }
0484   }
0485   for (int ilayer = 0; ilayer < 2; ilayer++)
0486   {
0487     for (int ieta = 0; ieta < _HCAL_NETA; ieta++)
0488     {
0489       for (int iphi = 0; iphi < _HCAL_NPHI; iphi++)
0490       {
0491         _TOWERMAP_STATUS_LAYER_ETA_PHI[ilayer][ieta][iphi] = -2;  // set tower does not exist
0492         _TOWERMAP_E_LAYER_ETA_PHI[ilayer][ieta][iphi] = 0;        // set zero energy
0493       }
0494     }
0495   }
0496 
0497   // setup
0498   std::vector<std::pair<int, float> > list_of_seeds;
0499 
0500   // translate towers to our internal representation
0501   if (_enable_EMCal)
0502   {
0503     TowerInfo *towerInfo = nullptr;
0504     unsigned int n_EM_towers = towerinfosEM->size();
0505     for (unsigned int iEM = 0; iEM < n_EM_towers; iEM++)
0506     {
0507       towerInfo = towerinfosEM->get_tower_at_channel(iEM);
0508       if (_only_good_towers && (!towerInfo->get_isGood()))
0509       {
0510         continue;
0511       }
0512       unsigned int towerinfo_key = towerinfosEM->encode_key(iEM);
0513       int ti_ieta = towerinfosEM->getTowerEtaBin(towerinfo_key);
0514       int ti_iphi = towerinfosEM->getTowerPhiBin(towerinfo_key);
0515       const RawTowerDefs::keytype key = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::CEMC, ti_ieta, ti_iphi);
0516 
0517       // RawTowerGeom *tower_geom = _geom_containers[2]->get_tower_geometry(key);
0518 
0519       // int ieta = _geom_containers[2]->get_etabin(tower_geom->get_eta());
0520       // int iphi = _geom_containers[2]->get_phibin(tower_geom->get_phi());
0521 
0522       int ieta = ti_ieta;
0523       int iphi = ti_iphi;
0524 
0525       float this_E = towerInfo->get_energy();
0526 
0527       // if not using abs E, short circuit all negative towers right here (same for IHCal, OHCal below)
0528       if (!_use_absE && this_E < 1.E-10)
0529       {
0530         continue;
0531       }
0532 
0533       _EMTOWERMAP_STATUS_ETA_PHI[ieta][iphi] = -1;  // change status to unknown
0534       _EMTOWERMAP_E_ETA_PHI[ieta][iphi] = this_E;
0535       _EMTOWERMAP_KEY_ETA_PHI[ieta][iphi] = key;
0536 
0537       // use fabs() here for simplicity - if we're not using abs E, negative towers are already excluded
0538       if (std::fabs(this_E) >= _sigma_seed * _noise_LAYER[2])
0539       {
0540         int ID = get_ID(2, ieta, iphi);
0541         list_of_seeds.emplace_back(ID, this_E);
0542         if (Verbosity() > 10)
0543         {
0544           std::cout << "RawClusterBuilderTopo::process_event: adding EMCal tower at ieta / iphi = " << ieta << " / " << iphi << " with E = " << this_E << std::endl;
0545           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;
0546         };
0547       }
0548     }
0549   }
0550 
0551   // translate towers to our internal representation
0552   if (_enable_HCal)
0553   {
0554     TowerInfo *towerInfo = nullptr;
0555     unsigned int n_IH_towers = towerinfosIH->size();
0556     for (unsigned int iIH = 0; iIH < n_IH_towers; iIH++)
0557     {
0558       towerInfo = towerinfosIH->get_tower_at_channel(iIH);
0559       if (_only_good_towers && (!towerInfo->get_isGood()))
0560       {
0561         continue;
0562       }
0563       unsigned int towerinfo_key = towerinfosIH->encode_key(iIH);
0564       int ti_ieta = towerinfosIH->getTowerEtaBin(towerinfo_key);
0565       int ti_iphi = towerinfosIH->getTowerPhiBin(towerinfo_key);
0566       const RawTowerDefs::keytype key = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::HCALIN, ti_ieta, ti_iphi);
0567 
0568       RawTowerGeom *tower_geom = _geom_containers[0]->get_tower_geometry(key);
0569 
0570       int ieta = _geom_containers[0]->get_etabin(tower_geom->get_eta());
0571       int iphi = _geom_containers[0]->get_phibin(tower_geom->get_phi());
0572       float this_E = towerInfo->get_energy();
0573 
0574       if (!_use_absE && this_E < 1.E-10)
0575       {
0576         continue;
0577       }
0578 
0579       _TOWERMAP_STATUS_LAYER_ETA_PHI[0][ieta][iphi] = -1;  // change status to unknown
0580       _TOWERMAP_E_LAYER_ETA_PHI[0][ieta][iphi] = this_E;
0581       _TOWERMAP_KEY_LAYER_ETA_PHI[0][ieta][iphi] = key;
0582 
0583       if (std::fabs(this_E) >= _sigma_seed * _noise_LAYER[0])
0584       {
0585         int ID = get_ID(0, ieta, iphi);
0586         list_of_seeds.emplace_back(ID, this_E);
0587         if (Verbosity() > 10)
0588         {
0589           std::cout << "RawClusterBuilderTopo::process_event: adding IHCal tower at ieta / iphi = " << ieta << " / " << iphi << " with E = " << this_E << std::endl;
0590           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;
0591         };
0592       }
0593     }
0594     unsigned int n_OH_towers = towerinfosOH->size();
0595     for (unsigned int iOH = 0; iOH < n_OH_towers; iOH++)
0596     {
0597       towerInfo = towerinfosOH->get_tower_at_channel(iOH);
0598       if (_only_good_towers && (!towerInfo->get_isGood()))
0599       {
0600         continue;
0601       }
0602       unsigned int towerinfo_key = towerinfosOH->encode_key(iOH);
0603       int ti_ieta = towerinfosOH->getTowerEtaBin(towerinfo_key);
0604       int ti_iphi = towerinfosOH->getTowerPhiBin(towerinfo_key);
0605       const RawTowerDefs::keytype key = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::HCALOUT, ti_ieta, ti_iphi);
0606 
0607       RawTowerGeom *tower_geom = _geom_containers[1]->get_tower_geometry(key);
0608 
0609       int ieta = _geom_containers[1]->get_etabin(tower_geom->get_eta());
0610       int iphi = _geom_containers[1]->get_phibin(tower_geom->get_phi());
0611       float this_E = towerInfo->get_energy();
0612 
0613       if (!_use_absE && this_E < 1.E-10)
0614       {
0615         continue;
0616       }
0617 
0618       _TOWERMAP_STATUS_LAYER_ETA_PHI[1][ieta][iphi] = -1;  // change status to unknown
0619       _TOWERMAP_E_LAYER_ETA_PHI[1][ieta][iphi] = this_E;
0620       _TOWERMAP_KEY_LAYER_ETA_PHI[1][ieta][iphi] = key;
0621 
0622       if (std::fabs(this_E) >= _sigma_seed * _noise_LAYER[1])
0623       {
0624         int ID = get_ID(1, ieta, iphi);
0625         list_of_seeds.emplace_back(ID, this_E);
0626         if (Verbosity() > 10)
0627         {
0628           std::cout << "RawClusterBuilderTopo::process_event: adding OHCal tower at ieta / iphi = " << ieta << " / " << iphi << " with E = " << this_E << std::endl;
0629           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;
0630         };
0631       }
0632     }
0633   }
0634 
0635   if (Verbosity() > 10)
0636   {
0637     for (unsigned int n = 0; n < list_of_seeds.size(); n++)
0638     {
0639       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;
0640     }
0641   }
0642 
0643   std::sort(list_of_seeds.begin(), list_of_seeds.end(), sort_by_pair_second);
0644 
0645   if (Verbosity() > 10)
0646   {
0647     for (unsigned int n = 0; n < list_of_seeds.size(); n++)
0648     {
0649       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;
0650     }
0651   }
0652 
0653   if (Verbosity() > 0)
0654   {
0655     std::cout << "RawClusterBuilderTopo::process_event: initialized with " << list_of_seeds.size() << " seeds with E > 4*sigma " << std::endl;
0656   }
0657 
0658   int cluster_index = 0;  // begin counting clusters
0659 
0660   std::vector<std::vector<int> > all_cluster_towers;  // store final cluster tower lists here
0661 
0662   while (!list_of_seeds.empty())
0663   {
0664     int seed_ID = list_of_seeds.at(0).first;
0665     list_of_seeds.erase(list_of_seeds.begin());
0666 
0667     if (Verbosity() > 5)
0668     {
0669       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;
0670     }
0671 
0672     // if this seed was already claimed by some other seed during its growth, remove it and do nothing
0673     int seed_status = get_status_from_ID(seed_ID);
0674     if (seed_status > -1)
0675     {
0676       if (Verbosity() > 10)
0677       {
0678         std::cout << " --> already owned by cluster # " << seed_status << std::endl;
0679       }
0680       continue;  // go onto the next iteration of the loop
0681     }
0682 
0683     // this seed tower now owned by new cluster
0684     set_status_by_ID(seed_ID, cluster_index);
0685 
0686     std::vector<int> cluster_tower_ID;
0687     cluster_tower_ID.push_back(seed_ID);
0688 
0689     std::vector<int> grow_tower_ID;
0690     grow_tower_ID.push_back(seed_ID);
0691 
0692     // iteratively process growth towers, adding > 2 * sigma neighbors to the list for further checking
0693 
0694     if (Verbosity() > 5)
0695     {
0696       std::cout << " RawClusterBuilderTopo::process_event: Entering Growth stage for cluster " << cluster_index << std::endl;
0697     }
0698 
0699     while (!grow_tower_ID.empty())
0700     {
0701       int grow_ID = grow_tower_ID.at(0);
0702       grow_tower_ID.erase(grow_tower_ID.begin());
0703 
0704       if (Verbosity() > 5)
0705       {
0706         std::cout << " --> cluster " << cluster_index << ", growth stage, examining neighbors of ID " << grow_ID << ", " << grow_tower_ID.size() << " grow towers left" << std::endl;
0707       }
0708 
0709       std::vector<int> adjacent_tower_IDs = get_adjacent_towers_by_ID(grow_ID);
0710 
0711       for (int this_adjacent_tower_ID : adjacent_tower_IDs)
0712       {
0713         if (Verbosity() > 10)
0714         {
0715           std::cout << " --> --> --> checking possible adjacent tower with ID " << this_adjacent_tower_ID << " : ";
0716         }
0717         int test_layer = get_ilayer_from_ID(this_adjacent_tower_ID);
0718 
0719         // if tower does not exist, continue
0720         if (get_status_from_ID(this_adjacent_tower_ID) == -2)
0721         {
0722           if (Verbosity() > 10)
0723           {
0724             std::cout << "does not exist " << std::endl;
0725           }
0726           continue;
0727         }
0728 
0729         // if tower is owned by THIS cluster already, continue
0730         if (get_status_from_ID(this_adjacent_tower_ID) == cluster_index)
0731         {
0732           if (Verbosity() > 10)
0733           {
0734             std::cout << "already owned by this cluster index " << cluster_index << std::endl;
0735           }
0736           continue;
0737         }
0738 
0739         // if tower has < 2*sigma energy, continue
0740         if (std::fabs(get_E_from_ID(this_adjacent_tower_ID)) < _sigma_grow * _noise_LAYER[test_layer])
0741         {
0742           if (Verbosity() > 10)
0743           {
0744             std::cout << "E = " << get_E_from_ID(this_adjacent_tower_ID) << " under 2*sigma threshold " << std::endl;
0745           }
0746           continue;
0747         }
0748 
0749         // if tower is owned by somebody else, continue (although should this really happen?)
0750         if (get_status_from_ID(this_adjacent_tower_ID) > -1)
0751         {
0752           if (Verbosity() > 10)
0753           {
0754             std::cout << "ERROR! in growth stage, encountered >2sigma tower which is already owned?!" << std::endl;
0755           }
0756           continue;
0757         }
0758 
0759         // tower good to be added to cluster and to list of grow towers
0760         grow_tower_ID.push_back(this_adjacent_tower_ID);
0761         cluster_tower_ID.push_back(this_adjacent_tower_ID);
0762         set_status_by_ID(this_adjacent_tower_ID, cluster_index);
0763         if (Verbosity() > 10)
0764         {
0765           std::cout << "add this tower ( ID " << this_adjacent_tower_ID << " ) to grow list " << std::endl;
0766         }
0767       }
0768 
0769       if (Verbosity() > 5)
0770       {
0771         std::cout << " --> after examining neighbors, grow list is now " << grow_tower_ID.size() << ", # of towers in cluster = " << cluster_tower_ID.size() << std::endl;
0772       }
0773     }
0774 
0775     // done growing cluster, now add on perimeter towers with E > 0 * sigma
0776     if (Verbosity() > 5)
0777     {
0778       std::cout << " RawClusterBuilderTopo::process_event: Entering Perimeter stage for cluster " << cluster_index << std::endl;
0779     }
0780     // we'll be adding on to the cluster list, so get the # of core towers first
0781     int n_core_towers = cluster_tower_ID.size();
0782 
0783     for (int ic = 0; ic < n_core_towers; ic++)
0784     {
0785       int core_ID = cluster_tower_ID.at(ic);
0786 
0787       if (Verbosity() > 5)
0788       {
0789         std::cout << " --> cluster " << cluster_index << ", perimeter stage, examining neighbors of ID " << core_ID << ", core cluster # " << ic << " of " << n_core_towers << " total " << std::endl;
0790       }
0791       std::vector<int> adjacent_tower_IDs = get_adjacent_towers_by_ID(core_ID);
0792 
0793       for (int this_adjacent_tower_ID : adjacent_tower_IDs)
0794       {
0795         if (Verbosity() > 10)
0796         {
0797           std::cout << " --> --> --> checking possible adjacent tower with ID " << this_adjacent_tower_ID << " : ";
0798         }
0799 
0800         int test_layer = get_ilayer_from_ID(this_adjacent_tower_ID);
0801 
0802         // if tower does not exist, continue
0803         if (get_status_from_ID(this_adjacent_tower_ID) == -2)
0804         {
0805           if (Verbosity() > 10)
0806           {
0807             std::cout << "does not exist " << std::endl;
0808           }
0809           continue;
0810         }
0811 
0812         // if tower is owned by somebody else (including current cluster), continue. ( allowed during perimeter fixing state )
0813         if (get_status_from_ID(this_adjacent_tower_ID) > -1)
0814         {
0815           if (Verbosity() > 10)
0816           {
0817             std::cout << "already owned by other cluster index " << get_status_from_ID(this_adjacent_tower_ID) << std::endl;
0818           }
0819           continue;
0820         }
0821 
0822         // if tower has < 0*sigma energy, continue
0823         if (std::fabs(get_E_from_ID(this_adjacent_tower_ID)) < _sigma_peri * _noise_LAYER[test_layer])
0824         {
0825           if (Verbosity() > 10)
0826           {
0827             std::cout << "E = " << get_E_from_ID(this_adjacent_tower_ID) << " under 0*sigma threshold " << std::endl;
0828           }
0829           continue;
0830         }
0831 
0832         // perimeter tower good to be added to cluster
0833         cluster_tower_ID.push_back(this_adjacent_tower_ID);
0834         set_status_by_ID(this_adjacent_tower_ID, cluster_index);
0835         if (Verbosity() > 10)
0836         {
0837           std::cout << "add this tower ( ID " << this_adjacent_tower_ID << " ) to cluster " << std::endl;
0838         }
0839       }
0840 
0841       if (Verbosity() > 5)
0842       {
0843         std::cout << " --> after examining perimeter neighbors, # of towers in cluster is now = " << cluster_tower_ID.size() << std::endl;
0844       }
0845     }
0846 
0847     // keep track of these
0848     all_cluster_towers.push_back(cluster_tower_ID);
0849 
0850     // increment cluster index for next one
0851     cluster_index++;
0852   }
0853 
0854   if (Verbosity() > 0)
0855   {
0856     std::cout << "RawClusterBuilderTopo::process_event: " << cluster_index << " topo-clusters initially reconstructed, entering splitting step" << std::endl;
0857   }
0858 
0859   int original_cluster_index = cluster_index;  // since it may be updated
0860 
0861   // now entering cluster splitting stage
0862 
0863   for (int cl = 0; cl < original_cluster_index; cl++)
0864   {
0865     std::vector<int> original_towers = all_cluster_towers.at(cl);
0866 
0867     if (!_do_split)
0868     {
0869       // don't run splitting, just export entire cluster as it is
0870       if (Verbosity() > 2)
0871       {
0872         std::cout << "RawClusterBuilderTopo::process_event: splitting step disabled, cluster " << cluster_index << " is final" << std::endl;
0873       }
0874       export_single_cluster(original_towers);
0875       continue;
0876     }
0877 
0878     std::vector<std::pair<int, float> > local_maxima_ID;
0879 
0880     // iterate through each tower, looking for maxima
0881     for (int tower_ID : original_towers)
0882     {
0883       if (Verbosity() > 10)
0884       {
0885         std::cout << " -> examining tower ID " << tower_ID << " for possible local maximum " << std::endl;
0886       }
0887 
0888       // check minimum energy
0889       if (get_E_from_ID(tower_ID) < _local_max_minE_LAYER[get_ilayer_from_ID(tower_ID)])
0890       {
0891         if (Verbosity() > 10)
0892         {
0893           std::cout << " -> -> energy E = " << get_E_from_ID(tower_ID) << " < " << _local_max_minE_LAYER[get_ilayer_from_ID(tower_ID)] << " too low" << std::endl;
0894         }
0895         continue;
0896       }
0897 
0898       // examine neighbors
0899       std::vector<int> adjacent_tower_IDs = get_adjacent_towers_by_ID(tower_ID);
0900       int neighbors_in_cluster = 0;
0901 
0902       // check for higher neighbor
0903       bool has_higher_neighbor = false;
0904       for (int this_adjacent_tower_ID : adjacent_tower_IDs)
0905       {
0906         if (get_status_from_ID(this_adjacent_tower_ID) != cl)
0907         {
0908           continue;  // only consider neighbors in cluster, obviously
0909         }
0910 
0911         neighbors_in_cluster++;
0912 
0913         if (get_E_from_ID(this_adjacent_tower_ID) > get_E_from_ID(tower_ID))
0914         {
0915           if (Verbosity() > 10)
0916           {
0917             std::cout << " -> -> has higher-energy neighbor ID / E = " << this_adjacent_tower_ID << " / " << get_E_from_ID(this_adjacent_tower_ID) << std::endl;
0918           }
0919           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
0920           break;
0921         }
0922       }
0923 
0924       if (has_higher_neighbor)
0925       {
0926         continue;  // if we broke out, now continue
0927       }
0928 
0929       // check number of neighbors
0930       if (neighbors_in_cluster < 4)
0931       {
0932         if (Verbosity() > 10)
0933         {
0934           std::cout << " -> -> too few neighbors N = " << neighbors_in_cluster << std::endl;
0935         }
0936         continue;
0937       }
0938 
0939       local_maxima_ID.emplace_back(tower_ID, get_E_from_ID(tower_ID));
0940     }
0941 
0942     // check for possible EMCal-OHCal seed overlaps
0943     for (unsigned int n = 0; n < local_maxima_ID.size(); n++)
0944     {
0945       // only look at I/OHCal local maxima
0946       std::pair<int, float> this_LM = local_maxima_ID.at(n);
0947       if (get_ilayer_from_ID(this_LM.first) == 2)
0948       {
0949         continue;
0950       }
0951 
0952       float this_phi = _geom_containers[get_ilayer_from_ID(this_LM.first)]->get_phicenter(get_iphi_from_ID(this_LM.first));
0953       if (this_phi > M_PI)
0954       {
0955         this_phi -= 2 * M_PI;
0956       }
0957       float this_eta = _geom_containers[get_ilayer_from_ID(this_LM.first)]->get_etacenter(get_ieta_from_ID(this_LM.first));
0958 
0959       bool has_EM_overlap = false;
0960 
0961       // check all other local maxima for overlaps
0962       for (unsigned int n2 = 0; n2 < local_maxima_ID.size(); n2++)
0963       {
0964         if (n == n2)
0965         {
0966           continue;  // don't check the same one
0967         }
0968 
0969         // only look at EMCal local mazima
0970         std::pair<int, float> this_LM2 = local_maxima_ID.at(n2);
0971         if (get_ilayer_from_ID(this_LM2.first) != 2)
0972         {
0973           continue;
0974         }
0975 
0976         float this_phi2 = _geom_containers[get_ilayer_from_ID(this_LM2.first)]->get_phicenter(get_iphi_from_ID(this_LM2.first));
0977         if (this_phi2 > M_PI)
0978         {
0979           this_phi -= 2 * M_PI;
0980         }
0981         float this_eta2 = _geom_containers[get_ilayer_from_ID(this_LM2.first)]->get_etacenter(get_ieta_from_ID(this_LM2.first));
0982 
0983         // calculate geometric dR
0984         float dR = calculate_dR(this_eta, this_eta2, this_phi, this_phi2);
0985 
0986         // check for and report overlaps
0987         if (dR < 0.15)
0988         {
0989           has_EM_overlap = true;
0990           if (Verbosity() > 2)
0991           {
0992             std::cout << "RawClusterBuilderTopo::process_event : removing I/OHal local maximum (ID,E,phi,eta = " << this_LM.first << ", " << this_LM.second << ", " << this_phi << ", " << this_eta << "), ";
0993             std::cout << "due to EM overlap (ID,E,phi,eta = " << this_LM2.first << ", " << this_LM2.second << ", " << this_phi2 << ", " << this_eta2 << "), dR = " << dR << std::endl;
0994           }
0995           break;
0996         }
0997       }
0998 
0999       if (has_EM_overlap)
1000       {
1001         // remove the I/OHCal local maximum from the list
1002         local_maxima_ID.erase(local_maxima_ID.begin() + n);
1003         // make sure to back up one index...
1004         n = n - 1;
1005       }  // otherwise, keep this local maximum
1006     }
1007 
1008     // only now print out full set of local maxima
1009     if (Verbosity() > 2)
1010     {
1011       for (auto this_LM : local_maxima_ID)
1012       {
1013         int tower_ID = this_LM.first;
1014         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) << ", ";
1015         float this_phi = _geom_containers[get_ilayer_from_ID(tower_ID)]->get_phicenter(get_iphi_from_ID(tower_ID));
1016         if (this_phi > M_PI)
1017         {
1018           this_phi -= 2 * M_PI;
1019         }
1020         std::cout << " eta / phi = " << _geom_containers[get_ilayer_from_ID(tower_ID)]->get_etacenter(get_ieta_from_ID(tower_ID)) << " / " << this_phi << std::endl;
1021       }
1022     }
1023 
1024     // do we have only 1 or 0 local maxima?
1025     if (local_maxima_ID.size() <= 1)
1026     {
1027       if (Verbosity() > 2)
1028       {
1029         std::cout << "RawClusterBuilderTopo::process_event cluster " << cl << " has only " << local_maxima_ID.size() << " local maxima, not splitting " << std::endl;
1030       }
1031       export_single_cluster(original_towers);
1032 
1033       continue;
1034     }
1035 
1036     // engage splitting procedure!
1037 
1038     if (Verbosity() > 2)
1039     {
1040       std::cout << "RawClusterBuilderTopo::process_event splitting cluster " << cl << " into " << local_maxima_ID.size() << " according to local maxima!" << std::endl;
1041     }
1042     // translate all cluster towers to a map which keeps track of their ownership
1043     // -1 means unseen
1044     // -2 means seen and in the seed list now (e.g. don't add it to the seed list again)
1045     // -3 shared tower, ignore going forward...
1046     std::map<int, std::pair<int, int> > tower_ownership;
1047     for (int &original_tower : original_towers)
1048     {
1049       tower_ownership[original_tower] = std::pair<int, int>(-1, -1);  // initialize all towers as un-seen
1050     }
1051     std::vector<int> seed_list;
1052     std::vector<int> neighbor_list;
1053     std::vector<int> shared_list;
1054 
1055     // sort maxima before populating seed list
1056     std::sort(local_maxima_ID.begin(), local_maxima_ID.end(), sort_by_pair_second);
1057 
1058     // initialize neighbor list
1059     for (unsigned int s = 0; s < local_maxima_ID.size(); s++)
1060     {
1061       tower_ownership[local_maxima_ID.at(s).first] = std::pair<int, int>(s, -1);
1062       neighbor_list.push_back(local_maxima_ID.at(s).first);
1063     }
1064 
1065     if (Verbosity() > 100)
1066     {
1067       for (int &original_tower : original_towers)
1068       {
1069         std::pair<int, int> the_pair = tower_ownership[original_tower];
1070         std::cout << " Debug Pre-Split: tower_ownership[ " << original_tower << " ] = ( " << the_pair.first << ", " << the_pair.second << " ) ";
1071         std::cout << " , layer / ieta / iphi = " << get_ilayer_from_ID(original_tower) << " / " << get_ieta_from_ID(original_tower) << " / " << get_iphi_from_ID(original_tower);
1072         std::cout << std::endl;
1073       }
1074     }
1075 
1076     bool first_pass = true;
1077 
1078     do
1079     {
1080       if (Verbosity() > 5)
1081       {
1082         std::cout << " -> starting split loop with " << seed_list.size() << " seed, " << neighbor_list.size() << " neighbor, and " << shared_list.size() << " shared towers " << std::endl;
1083       }
1084       // go through neighbor list, assigning ownership only via the seed list
1085       std::vector<int> new_ownerships;
1086 
1087       for (unsigned int n = 0; n < neighbor_list.size(); n++)
1088       {
1089         int neighbor_ID = neighbor_list.at(n);
1090 
1091         if (Verbosity() > 10)
1092         {
1093           std::cout << " -> -> looking at neighbor " << n << " (tower ID " << neighbor_ID << " ) of " << neighbor_list.size() << " total" << std::endl;
1094         }
1095         if (first_pass)
1096         {
1097           if (Verbosity() > 10)
1098           {
1099             std::cout << " -> -> -> special first pass rules, this tower already owned by pseudocluster " << tower_ownership[neighbor_ID].first << std::endl;
1100           }
1101           new_ownerships.push_back(tower_ownership[neighbor_ID].first);
1102         }
1103         else
1104         {
1105           std::map<int, bool> pseudocluster_adjacency;
1106           for (unsigned int s = 0; s < local_maxima_ID.size(); s++)
1107           {
1108             pseudocluster_adjacency[s] = false;
1109           }
1110           // look over all towers THIS one is adjacent to, and count up...
1111           std::vector<int> adjacent_tower_IDs = get_adjacent_towers_by_ID(neighbor_ID);
1112 
1113           for (int this_adjacent_tower_ID : adjacent_tower_IDs)
1114           {
1115             if (get_status_from_ID(this_adjacent_tower_ID) != cl)
1116             {
1117               continue;
1118             }
1119 
1120             if (tower_ownership[this_adjacent_tower_ID].first > -1)
1121             {
1122               if (Verbosity() > 20)
1123               {
1124                 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;
1125               }
1126               pseudocluster_adjacency[tower_ownership[this_adjacent_tower_ID].first] = true;
1127             }
1128           }
1129           int n_pseudocluster_adjacent = 0;
1130           int last_adjacent_pseudocluster = -1;
1131           for (unsigned int s = 0; s < local_maxima_ID.size(); s++)
1132           {
1133             if (pseudocluster_adjacency[s])
1134             {
1135               last_adjacent_pseudocluster = s;
1136               n_pseudocluster_adjacent++;
1137               if (Verbosity() > 20)
1138               {
1139                 std::cout << " -> -> adjacent to pseudocluster " << s << std::endl;
1140               }
1141             }
1142           }
1143 
1144           if (n_pseudocluster_adjacent == 0)
1145           {
1146             std::cout << " -> -> ERROR! How can a neighbor tower at this stage be adjacent to no pseudoclusters?? " << std::endl;
1147             new_ownerships.push_back(9999);
1148           }
1149           else if (n_pseudocluster_adjacent == 1)
1150           {
1151             if (Verbosity() > 10)
1152             {
1153               std::cout << " -> -> neighbor tower " << neighbor_ID << " is ONLY adjacent to one pseudocluster # " << last_adjacent_pseudocluster << std::endl;
1154             }
1155             new_ownerships.push_back(last_adjacent_pseudocluster);
1156           }
1157           else
1158           {
1159             if (Verbosity() > 10)
1160             {
1161               std::cout << " -> -> neighbor tower " << neighbor_ID << " is adjacent to " << n_pseudocluster_adjacent << " pseudoclusters, move to shared list " << std::endl;
1162             }
1163             new_ownerships.push_back(-3);
1164           }
1165         }
1166       }
1167 
1168       if (Verbosity() > 5)
1169       {
1170         std::cout << " -> now updating status of all " << neighbor_list.size() << " original neighbors " << std::endl;
1171       }
1172       // transfer neighbor list to seed list or shared list
1173       for (unsigned int n = 0; n < neighbor_list.size(); n++)
1174       {
1175         int neighbor_ID = neighbor_list.at(n);
1176         if (new_ownerships.at(n) > -1)
1177         {
1178           tower_ownership[neighbor_ID] = std::pair<int, int>(new_ownerships.at(n), -1);
1179           seed_list.push_back(neighbor_ID);
1180           if (Verbosity() > 20)
1181           {
1182             std::cout << " -> -> neighbor ID " << neighbor_ID << " has new status " << new_ownerships.at(n) << std::endl;
1183           }
1184         }
1185         if (new_ownerships.at(n) == -3)
1186         {
1187           tower_ownership[neighbor_ID] = std::pair<int, int>(-3, -1);
1188           shared_list.push_back(neighbor_ID);
1189           if (Verbosity() > 20)
1190           {
1191             std::cout << " -> -> neighbor ID " << neighbor_ID << " has new status " << -3 << std::endl;
1192           }
1193         }
1194       }
1195 
1196       if (Verbosity() > 5)
1197       {
1198         std::cout << " producing a new neighbor list ... " << std::endl;
1199       }
1200       // populate a new neighbor list from the about-to-be-owned towers before transferring this one
1201       std::list<int> new_neighbor_list;
1202       for (unsigned int n = 0; n < neighbor_list.size(); n++)
1203       {
1204         int neighbor_ID = neighbor_list.at(n);
1205         if (new_ownerships.at(n) > -1)
1206         {
1207           std::vector<int> adjacent_tower_IDs = get_adjacent_towers_by_ID(neighbor_ID);
1208 
1209           for (int this_adjacent_tower_ID : adjacent_tower_IDs)
1210           {
1211             if (get_status_from_ID(this_adjacent_tower_ID) != cl)
1212             {
1213               continue;
1214             }
1215             if (tower_ownership[this_adjacent_tower_ID].first == -1)
1216             {
1217               new_neighbor_list.push_back(this_adjacent_tower_ID);
1218               if (Verbosity() > 5)
1219               {
1220                 std::cout << " -> queueing up to add tower " << this_adjacent_tower_ID << " , neighbor of tower " << neighbor_ID << " to new neighbor list" << std::endl;
1221               }
1222             }
1223           }
1224         }
1225       }
1226 
1227       if (Verbosity() > 5)
1228       {
1229         std::cout << " new neighbor list has size " << new_neighbor_list.size() << ", but after removing duplicate elements: ";
1230       }
1231 
1232       new_neighbor_list.sort();
1233       new_neighbor_list.unique();
1234 
1235       if (Verbosity() > 5)
1236       {
1237         std::cout << new_neighbor_list.size() << std::endl;
1238       }
1239 
1240       // clear neighbor list!
1241       neighbor_list.clear();
1242 
1243       // now transfer over new neighbor list
1244       for (; !new_neighbor_list.empty();)
1245       {
1246         neighbor_list.push_back(new_neighbor_list.front());
1247         new_neighbor_list.pop_front();
1248       }
1249 
1250       first_pass = false;
1251 
1252     } while (!neighbor_list.empty());
1253 
1254     if (Verbosity() > 100)
1255     {
1256       for (int &original_tower : original_towers)
1257       {
1258         std::pair<int, int> the_pair = tower_ownership[original_tower];
1259         std::cout << " Debug Mid-Split: tower_ownership[ " << original_tower << " ] = ( " << the_pair.first << ", " << the_pair.second << " ) ";
1260         std::cout << " , layer / ieta / iphi = " << get_ilayer_from_ID(original_tower) << " / " << get_ieta_from_ID(original_tower) << " / " << get_iphi_from_ID(original_tower);
1261         std::cout << std::endl;
1262         if (the_pair.first == -1)
1263         {
1264           std::vector<int> adjacent_tower_IDs = get_adjacent_towers_by_ID(original_tower);
1265 
1266           for (int this_adjacent_tower_ID : adjacent_tower_IDs)
1267           {
1268             if (get_status_from_ID(this_adjacent_tower_ID) != cl)
1269             {
1270               continue;
1271             }
1272             std::cout << "    -> adjacent to add tower " << this_adjacent_tower_ID << " , which has status " << tower_ownership[this_adjacent_tower_ID].first << std::endl;
1273           }
1274         }
1275       }
1276     }
1277 
1278     // calculate pseudocluster energies and positions
1279     std::vector<float> pseudocluster_sumeta;
1280     std::vector<float> pseudocluster_sumphi;
1281     std::vector<float> pseudocluster_sumE;
1282     std::vector<int> pseudocluster_ntower;
1283     std::vector<float> pseudocluster_eta;
1284     std::vector<float> pseudocluster_phi;
1285 
1286     pseudocluster_sumeta.resize(local_maxima_ID.size(), 0);
1287     pseudocluster_sumphi.resize(local_maxima_ID.size(), 0);
1288     pseudocluster_sumE.resize(local_maxima_ID.size(), 0);
1289     pseudocluster_ntower.resize(local_maxima_ID.size(), 0);
1290 
1291     for (int &original_tower : original_towers)
1292     {
1293       std::pair<int, int> the_pair = tower_ownership[original_tower];
1294       if (the_pair.first > -1)
1295       {
1296         float this_ID = original_tower;
1297         pseudocluster_sumE[the_pair.first] += get_E_from_ID(this_ID);
1298         float this_eta = _geom_containers[get_ilayer_from_ID(this_ID)]->get_etacenter(get_ieta_from_ID(this_ID));
1299         float this_phi = _geom_containers[get_ilayer_from_ID(this_ID)]->get_phicenter(get_iphi_from_ID(this_ID));
1300 
1301         pseudocluster_sumeta[the_pair.first] += this_eta;
1302         pseudocluster_sumphi[the_pair.first] += this_phi;
1303         pseudocluster_ntower[the_pair.first] += 1;
1304       }
1305     }
1306 
1307     for (unsigned int pc = 0; pc < local_maxima_ID.size(); pc++)
1308     {
1309       pseudocluster_eta.push_back(pseudocluster_sumeta.at(pc) / pseudocluster_ntower.at(pc));
1310       pseudocluster_phi.push_back(pseudocluster_sumphi.at(pc) / pseudocluster_ntower.at(pc));
1311 
1312       if (Verbosity() > 2)
1313       {
1314         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;
1315       }
1316     }
1317 
1318     if (Verbosity() > 2)
1319     {
1320       std::cout << "RawClusterBuilderTopo::process_event now splitting up shared clusters (including unassigned clusters), initial shared list has size " << shared_list.size() << std::endl;
1321     }
1322     // iterate through shared cells, identifying which two they belong to
1323     while (!shared_list.empty())
1324     {
1325       // pick the first cell and pop off list
1326       int shared_ID = shared_list.at(0);
1327       shared_list.erase(shared_list.begin());
1328 
1329       if (Verbosity() > 5)
1330       {
1331         std::cout << " -> looking at shared tower " << shared_ID << ", after this one there are " << shared_list.size() << " shared towers left " << std::endl;
1332       }
1333       // look through adjacent pseudoclusters, taking two with highest energies
1334       std::vector<bool> pseudocluster_adjacency;
1335       pseudocluster_adjacency.resize(local_maxima_ID.size(), false);
1336 
1337       std::vector<int> adjacent_tower_IDs = get_adjacent_towers_by_ID(shared_ID);
1338 
1339       for (int this_adjacent_tower_ID : adjacent_tower_IDs)
1340       {
1341         if (get_status_from_ID(this_adjacent_tower_ID) != cl)
1342         {
1343           continue;
1344         }
1345         if (tower_ownership[this_adjacent_tower_ID].first > -1)
1346         {
1347           pseudocluster_adjacency[tower_ownership[this_adjacent_tower_ID].first] = true;
1348         }
1349         if (tower_ownership[this_adjacent_tower_ID].second > -1)
1350         {  // can inherit adjacency from shared cluster
1351           pseudocluster_adjacency[tower_ownership[this_adjacent_tower_ID].second] = true;
1352         }
1353         // at the same time, add unowned towers to the list for later examination
1354         if (tower_ownership[this_adjacent_tower_ID].first == -1)
1355         {
1356           shared_list.push_back(this_adjacent_tower_ID);
1357           tower_ownership[this_adjacent_tower_ID] = std::pair<int, int>(-3, -1);
1358           if (Verbosity() > 10)
1359           {
1360             std::cout << " -> while looking at neighbors, have added un-examined tower " << this_adjacent_tower_ID << " to shared list " << std::endl;
1361           }
1362         }
1363       }
1364 
1365       // now figure out which pseudoclusters this shared tower is adjacent to...
1366       int highest_pseudocluster_index = -1;
1367       int second_highest_pseudocluster_index = -1;
1368 
1369       float highest_pseudocluster_E = -999;
1370       float second_highest_pseudocluster_E = -999;
1371 
1372       for (unsigned int n = 0; n < pseudocluster_adjacency.size(); n++)
1373       {
1374         if (!pseudocluster_adjacency[n])
1375         {
1376           continue;
1377         }
1378 
1379         if (pseudocluster_sumE[n] > highest_pseudocluster_E)
1380         {
1381           second_highest_pseudocluster_E = highest_pseudocluster_E;
1382           second_highest_pseudocluster_index = highest_pseudocluster_index;
1383 
1384           highest_pseudocluster_E = pseudocluster_sumE[n];
1385           highest_pseudocluster_index = n;
1386         }
1387         else if (pseudocluster_sumE[n] > second_highest_pseudocluster_E)
1388         {
1389           second_highest_pseudocluster_E = pseudocluster_sumE[n];
1390           second_highest_pseudocluster_index = n;
1391         }
1392       }
1393 
1394       if (Verbosity() > 5)
1395       {
1396         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;
1397       }
1398       // assign these clusters as owners
1399       tower_ownership[shared_ID] = std::pair<int, int>(highest_pseudocluster_index, second_highest_pseudocluster_index);
1400     }
1401 
1402     if (Verbosity() > 100)
1403     {
1404       for (int &original_tower : original_towers)
1405       {
1406         std::pair<int, int> the_pair = tower_ownership[original_tower];
1407         std::cout << " Debug Post-Split: tower_ownership[ " << original_tower << " ] = ( " << the_pair.first << ", " << the_pair.second << " ) ";
1408         std::cout << " , layer / ieta / iphi = " << get_ilayer_from_ID(original_tower) << " / " << get_ieta_from_ID(original_tower) << " / " << get_iphi_from_ID(original_tower);
1409         std::cout << std::endl;
1410         if (the_pair.first == -1)
1411         {
1412           std::vector<int> adjacent_tower_IDs = get_adjacent_towers_by_ID(original_tower);
1413 
1414           for (int this_adjacent_tower_ID : adjacent_tower_IDs)
1415           {
1416             if (get_status_from_ID(this_adjacent_tower_ID) != cl)
1417             {
1418               continue;
1419             }
1420             std::cout << " -> adjacent to add tower " << this_adjacent_tower_ID << " , which has status " << tower_ownership[this_adjacent_tower_ID].first << std::endl;
1421           }
1422         }
1423       }
1424     }
1425 
1426     // call helper function
1427     export_clusters(original_towers, tower_ownership, local_maxima_ID.size(), pseudocluster_sumE, pseudocluster_eta, pseudocluster_phi);
1428   }
1429 
1430   if (Verbosity() > 1)
1431   {
1432     std::cout << "RawClusterBuilderTopo::process_event after splitting (if any) final clusters output to node are: " << std::endl;
1433     RawClusterContainer::ConstRange begin_end = _clusters->getClusters();
1434     int ncl = 0;
1435     for (RawClusterContainer::ConstIterator hiter = begin_end.first; hiter != begin_end.second; ++hiter)
1436     {
1437       std::cout << "-> #" << ncl++ << " ";
1438       hiter->second->identify();
1439       std::cout << std::endl;
1440     }
1441   }
1442 
1443   return Fun4AllReturnCodes::EVENT_OK;
1444 }
1445 
1446 int RawClusterBuilderTopo::End(PHCompositeNode * /*topNode*/)
1447 {
1448   return Fun4AllReturnCodes::EVENT_OK;
1449 }
1450 
1451 void RawClusterBuilderTopo::CreateNodes(PHCompositeNode *topNode)
1452 {
1453   PHNodeIterator iter(topNode);
1454 
1455   PHCompositeNode *dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
1456   if (!dstNode)
1457   {
1458     std::cerr << PHWHERE << "DST Node missing, doing nothing." << std::endl;
1459     throw std::runtime_error("Failed to find DST node in RawClusterBuilderTopo::CreateNodes");
1460   }
1461 
1462   PHNodeIterator dstiter(dstNode);
1463   PHCompositeNode *DetNode = dynamic_cast<PHCompositeNode *>(dstiter.findFirst("PHCompositeNode", "TOPOCLUSTER"));
1464   if (!DetNode)
1465   {
1466     DetNode = new PHCompositeNode("TOPOCLUSTER");
1467     dstNode->addNode(DetNode);
1468   }
1469 
1470   _clusters = new RawClusterContainer();
1471 
1472   PHIODataNode<PHObject> *clusterNode = new PHIODataNode<PHObject>(_clusters, ClusterNodeName, "PHObject");
1473   DetNode->addNode(clusterNode);
1474 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!