trigger::IdealClusterBuilder Class Reference

Public Member Functions
void	AddHit (Hit h)
std::vector< Cluster >	GetClusters ()
void	SetClusterGeo (ClusterGeometry *_g)
virtual void	BuildClusters ()
std::vector< Cluster >	Build2dClustersLayer (std::vector< Hit > hits)
void	Build2dClusters ()
void	Build3dClusters ()
void	Fit (Cluster &c3)

Public Attributes
std::vector< Hit >	all_hits {}
std::vector< Cluster >	all_clusters {}
ClusterGeometry *	g
float	seed_thresh = 0
float	neighb_thresh = 0
int	n_neighbors = 1
bool	split_energy = true
bool	use_towers = false
const int	LAYER_MAX = 35
const int	LAYER_SHOWERMAX = 7
const int	LAYER_SEEDMIN = 3
const int	LAYER_SEEDMAX = 15
const float	MIN_TP_ENERGY = 0.5
int	DEPTH_GOOD = 5
bool	debug = false

Detailed Description

Definition at line 142 of file IdealClusterBuilder.h.

Member Function Documentation

AddHit()

void trigger::IdealClusterBuilder::AddHit ( Hit  h )

inline

Definition at line 166 of file IdealClusterBuilder.h.

                     {
    if (h.layer >= LAYER_MAX) return;
    auto p = std::make_pair(h.cell_id, h.module_id);
    h.id = g->reverse_id_map[p];
    all_hits.push_back(h);
  }

Build2dClusters()

void trigger::IdealClusterBuilder::Build2dClusters

(

)

Definition at line 207 of file IdealClusterBuilder.cxx.

                                          {
  // first partition hits by layer
  std::map<int, std::vector<Hit> > layer_hits;  // id(xy) to Hit
  for (const auto hit : all_hits) {
    auto l = hit.layer;
    if (layer_hits.count(l)) {
      layer_hits[l].push_back(hit);
    } else {
      layer_hits[l] = {hit};
    }
  }
 
  // run clustering in each layer and add to the list
  for (auto &pair : layer_hits) {
    if (debug) {
      cout << "Found " << pair.second.size() << " hits in layer " << pair.first
           << endl;
    }
    auto clus = Build2dClustersLayer(pair.second);
    all_clusters.insert(all_clusters.end(), clus.begin(), clus.end());
  }
}

Build2dClustersLayer()

std::vector< Cluster > trigger::IdealClusterBuilder::Build2dClustersLayer

(

std::vector< Hit >

hits

)

Definition at line 58 of file IdealClusterBuilder.cxx.

                         {
  // Re-index by id
  std::map<int, Hit> hits_by_id;
  for (auto &hit : hits) hits_by_id[hit.id] = hit;
 
  if (debug) {
    cout << "--------\nBuild2dClustersLayer Input Hits" << endl;
    for (auto &hitpair : hits_by_id) hitpair.second.Print();
  }
 
  // Find seeds
  std::vector<Cluster> clusters;
  for (auto &hitpair : hits_by_id) {
    auto &hit = hitpair.second;
    bool isLocalMax = true;
    for (auto n : g->neighbors[hit.id]) {
      // cout << "  checking " << n << endl;
      if (hits_by_id.count(n) && hits_by_id[n].e > hit.e) isLocalMax = false;
    }
    // if(debug) cout << hit.e << " " << hit.id << " "
    //  << hit.layer << " isMax=" << isLocalMax << endl;
    if (isLocalMax && (hit.e > seed_thresh)) {
      hit.used = true;
      Cluster c;
      c.hits.push_back(hit);
      c.e = hit.e;
      c.x = hit.x;
      c.y = hit.y;
      c.z = hit.z;
      c.seed = hit.id;
      c.module = g->id_map[hit.id].second;
      c.layer = hit.layer;
      clusters.push_back(c);
    }
  }
 
  if (debug) {
    cout << "--------\nAfter seed-finding" << endl;
    for (auto &hitpair : hits_by_id) hitpair.second.Print();
    for (auto &c : clusters) c.Print();
  }
 
  // Add neighbors up to the specified limit
  int i_neighbor = 0;
  while (i_neighbor < n_neighbors) {
    // find (unused) neighbors for all clusters
    std::map<int, std::vector<int> > assoc_clus2hitIDs;
    for (int iclus = 0; iclus < clusters.size(); iclus++) {
      auto &clus = clusters[iclus];
      std::vector<int> neighbors;
      for (const auto &hit : clus.hits) {
        for (auto n : g->neighbors[hit.id]) {
          if (hits_by_id.count(n) && !hits_by_id[n].used &&
              hits_by_id[n].e > neighb_thresh) {
            neighbors.push_back(n);
          }
        }
      }
      assoc_clus2hitIDs[iclus] = neighbors;
    }
 
    // check how many clusters to which each hit is assoc
    std::map<int, std::vector<int> > assoc_hitID2clusters;
    for (auto clus2hitID : assoc_clus2hitIDs) {
      auto iclus = clus2hitID.first;
      auto &hitIDs = clus2hitID.second;
      for (const auto &hitID : hitIDs) {
        if (assoc_hitID2clusters.count(hitID))
          assoc_hitID2clusters[hitID].push_back(iclus);
        else
          assoc_hitID2clusters[hitID] = {iclus};
      }
    }
 
    // add associated hits to clusters
    //   (w/ optional e-splitting)
    for (auto hitID2clusters : assoc_hitID2clusters) {
      auto hitID = hitID2clusters.first;
      auto iclusters = hitID2clusters.second;
      if (iclusters.size() == 1) {
        // simply add cell to the cluster
        auto &hit = hits_by_id[hitID];
        auto iclus = iclusters[0];
        hit.used = true;
        clusters[iclus].hits.push_back(hit);
        clusters[iclus].e += hit.e;
      } else {
        auto &hit = hits_by_id[hitID];
        hit.used = true;
        float esum = 0;
        for (auto iclus : iclusters) {
          esum += clusters[iclus].e;
        }
        for (auto iclus : iclusters) {
          Hit newHit = hit;
          if (split_energy) newHit.e = hit.e * clusters[iclus].e / esum;
          clusters[iclus].hits.push_back(newHit);
          clusters[iclus].e += newHit.e;
        }
      }
    }
 
    // rebuild the clusters and return
    for (auto &c : clusters) {
      c.e = 0;
      c.x = 0;
      c.y = 0;
      c.z = 0;
      c.xx = 0;
      c.yy = 0;
      c.zz = 0;
      float sumw = 0;
      for (auto hit : c.hits) {
        // if(debug) cout << hit.x << " " << hit.y << " " << hit.z << endl;
        c.e += hit.e;
        // cout << "2d: " << h.e << " " << log(h.e/MIN_TP_ENERGY) << endl;
        float w = std::max(0., log(hit.e / MIN_TP_ENERGY));  // use log-e wgt
        c.x += hit.x * w;
        c.y += hit.y * w;
        c.z += hit.z * w;
        c.xx += hit.x * hit.x * w;
        c.yy += hit.y * hit.y * w;
        c.zz += hit.z * hit.z * w;
        sumw += w;
      }
      c.x /= sumw;
      c.y /= sumw;
      c.z /= sumw;
      c.xx /= sumw;  // now is <x^2>
      c.yy /= sumw;
      c.zz /= sumw;
      c.xx = sqrt(c.xx - c.x * c.x);  // now is sqrt(<x^2>-<x>^2)
      c.yy = sqrt(c.yy - c.y * c.y);
      c.zz = sqrt(c.zz - c.z * c.z);
    }
 
    i_neighbor++;
 
    if (debug) {
      cout << "--------\nAfter " << i_neighbor << " neighbors" << endl;
      for (auto &hitpair : hits_by_id) hitpair.second.Print();
      for (auto &c : clusters) c.Print();
    }
  }
 
  return clusters;
}

Build3dClusters()

void trigger::IdealClusterBuilder::Build3dClusters

(

)

Definition at line 230 of file IdealClusterBuilder.cxx.

                                          {
  if (debug) {
    cout << "--------\nBuilding 3d clusters" << endl;
  }
 
  // first partition 2d clusters by layer
  std::vector<std::vector<Cluster> > layer_clusters;
  layer_clusters.resize(LAYER_MAX);  // first 20 layers
  for (auto &clus : all_clusters) {
    layer_clusters[clus.layer].push_back(clus);
  }
 
  // sort by layer
  for (auto &clusters : layer_clusters) ESort(clusters);
 
  if (debug) {
    cout << "--------\n3d: sorted 2d inputs" << endl;
    for (auto &clusters : layer_clusters)
      for (auto &c : clusters) c.Print(g);
  }
 
  // Pass through clusters from layer 0 to last,
  //   starting with highest energy
  bool building = true;
  std::vector<Cluster> clusters3d;
  while (building) {
    // find the seed cluster
    Cluster cluster3d;
    cluster3d.is2D = false;
    cluster3d.first_layer = LAYER_SHOWERMAX;
    cluster3d.last_layer = LAYER_SHOWERMAX;
    int test_layer;
    for (int ilayer = 0; ilayer < LAYER_MAX; ilayer++) {
      if (LAYER_SHOWERMAX + ilayer < LAYER_MAX) {
        // walk to back of ECal from shower max
        test_layer = LAYER_SHOWERMAX + ilayer;  // 7,8,9,...19
      } else {
        // then to front of ECal
        test_layer = LAYER_MAX - ilayer - 1;  // 20-13-1=6,5,4,...
      }
 
      auto &clusters2d = layer_clusters[test_layer];
 
      // still must find the 3d seed
      if (cluster3d.depth == 0) {
        // only seed from the middle layers
        if (test_layer > LAYER_SEEDMAX || test_layer < LAYER_SEEDMIN) continue;
        if (clusters2d.size()) {
          if (debug) {
            cout << " 3d seed: ";
            clusters2d[0].Print(g);
          }
          // cout << "got seed" << endl;
          //  found seed
          cluster3d.clusters2d = {clusters2d[0]};
          cluster3d.first_layer = test_layer;
          cluster3d.last_layer = test_layer;
          cluster3d.depth = 1;
          // remove (first) 2d cluster from list
          clusters2d.erase(clusters2d.begin());
        }
      } else {
        // looking to extend the 3d seed
        // grow if 2d seed is a neighbor
        // auto &last_seed2d = clusters3d.back().seed;
        auto &last_seed2d = cluster3d.clusters2d.back().seed;
        // case where we begin extending cluster backward->forward
        if (test_layer == LAYER_SHOWERMAX - 1)
          last_seed2d = cluster3d.clusters2d.front().seed;
        if (debug) {
          //    cout << cluster3d.clusters2d.size() << endl;
          //    cluster3d.clusters2d.back().Print();
          cout << "  check 3d w/ seed id " << last_seed2d << endl;
          // cout << "   from #cands: " << clusters2d.size() << endl;
        }
        for (int iclus2d = 0; iclus2d < clusters2d.size(); iclus2d++) {
          if (debug) {
            cout << "  -- " << iclus2d << endl;
            //     clusters2d[iclus2d].Print(g);
            // cout << "  check ext " << seed2d << endl;
          }
          auto &seed2d = clusters2d[iclus2d].seed;
          if (last_seed2d == seed2d || g->CheckNeighbor(last_seed2d, seed2d)) {
            // if(debug){
            //  cout << " extend: ";
            //  clusters2d[iclus2d].Print(g);
            // }
            // add to 3d cluster
            cluster3d.clusters2d.push_back(clusters2d[iclus2d]);
            cluster3d.depth++;
            if (test_layer < cluster3d.first_layer)
              cluster3d.first_layer = test_layer;
            if (test_layer > cluster3d.last_layer)
              cluster3d.last_layer = test_layer;
            // remove from list
            clusters2d.erase(clusters2d.begin() + iclus2d);
            // proceed to next layer
            break;
          }
        }
      }
    }
    // done with all layers. finish or store cluster
    if (cluster3d.depth == 0) {
      building = false;
    } else {
      // cout << "storing 3d cluster" << endl;
      if (cluster3d.depth >= DEPTH_GOOD) clusters3d.push_back(cluster3d);
    }
  }
 
  // post-process 3d clusters here
  for (auto &c : clusters3d) {
    c.e = 0;
    c.x = 0;
    c.y = 0;
    c.z = 0;
    c.xx = 0;
    c.yy = 0;
    c.zz = 0;
    float sumw = 0;
    for (auto &c2 : c.clusters2d) {
      c.e += c2.e;
      // cout << "3d: " << c2.e << " " << log(c2.e/MIN_TP_ENERGY) << endl;
      float w = std::max(0., log(c2.e / MIN_TP_ENERGY));  // use log-e wgt
      c.x += c2.x * w;
      c.y += c2.y * w;
      c.z += c2.z * w;
      c.xx += c2.x * c2.x * w;
      c.yy += c2.y * c2.y * w;
      c.zz += c2.z * c2.z * w;
      sumw += w;
    }
    // cout << "sum: " << sumw << endl;
    // cout << "x: " << c.x << endl;
    c.x /= sumw;
    // cout << "x: " << c.x << endl;
    c.y /= sumw;
    c.z /= sumw;
    c.xx /= sumw;  // now is <x^2>
    c.yy /= sumw;
    c.zz /= sumw;
    c.xx = sqrt(c.xx - c.x * c.x);  // now is sqrt(<x^2>-<x>^2)
    c.yy = sqrt(c.yy - c.y * c.y);
    c.zz = sqrt(c.zz - c.z * c.z);
    Fit(c);  // calc dx/dz, dy/dz
  }
 
  if (debug) {
    cout << "--------\nFound 3d clusters" << endl;
    for (auto &c : clusters3d) c.Print3d();
  }
 
  // std::map<int, std::vector<Cluster> > layer_clusters; // id(xy) to Hit
  // for(const auto clus : all_clusters){
  //     auto l = clus.layer;
  //     if (layer_clusters.count(l)){
  //    layer_clusters[l].push_back(clus);
  //     } else {
  //    layer_clusters[l]={clus};
  //     }
  // }
  // // sort by layer
  // for(auto &pair : layer_clusters){
  //     auto &clusters = pair.second;
 
  //     if(debug && clusters.size()>2){
  //    cout << "--------\nBefore sort " << endl;
  //    for(auto &c : clusters) c.Print();
  //     }
  //     ESort(clusters);
  //     if(debug && clusters.size()>2){
  //    cout << "-- After sort " << endl;
  //    for(auto &c : clusters) c.Print();
  //     }
  // }
 
  all_clusters.clear();
  all_clusters.insert(all_clusters.begin(), clusters3d.begin(),
                      clusters3d.end());
}

BuildClusters()

void trigger::IdealClusterBuilder::BuildClusters ( )

virtual

Definition at line 412 of file IdealClusterBuilder.cxx.

                                        {
  if (debug) {
    cout << "--------\nAll hits" << endl;
    for (auto &hit : all_hits) hit.Print();
  }
 
  if (use_towers) {
    // project hits in z to form towers
    std::map<int, Hit> towers;  // id(xy) to Hit
    for (const auto hit : all_hits) {
      if (towers.count(hit.id)) {
        towers[hit.id].e += hit.e;
        towers[hit.id].nSubHit++;
      } else {
        towers[hit.id] = hit;
        towers[hit.id].layer = 0;
        towers[hit.id].z = 0;
        towers[hit.id].nSubHit = 1;
      }
    }
    all_hits.clear();
    for (const auto t : towers) all_hits.push_back(t.second);
 
    if (debug) {
      cout << "--------\nHits after towers" << endl;
      for (auto &hit : all_hits) hit.Print();
    }
  }
 
  // Cluster the hits in each plane
  Build2dClusters();
 
  if (!use_towers) {
    Build3dClusters();
  }
  ESort(all_clusters);
};

Fit()

void trigger::IdealClusterBuilder::Fit

(

Cluster &

c3

)

Definition at line 450 of file IdealClusterBuilder.cxx.

                                         {
  // TODO: think about whether to incorporate uncertainties
  //   into the fit (RMSs), or weight each layer in the fit.
 
  // skip short clusters
  if (c3.clusters2d.size() < 4) return;
 
  // std::vector logE;
  std::vector<float> x;
  std::vector<float> y;
  std::vector<float> z;
  for (const auto &c2 : c3.clusters2d) {
    // logE.push_back( log(c2.e) );
    x.push_back(c2.x);
    y.push_back(c2.y);
    z.push_back(c2.z);
  }
  TGraph gxz(z.size(), z.data(), x.data());
  auto r_xz = gxz.Fit("pol1", "SQ");  // p0 + x*p1
  c3.dxdz = r_xz->Value(1);
  c3.dxdze = r_xz->ParError(1);
 
  TGraph gyz(z.size(), z.data(), y.data());
  auto r_yz = gyz.Fit("pol1", "SQ");  // p0 + x*p1
  c3.dydz = r_yz->Value(1);
  c3.dydze = r_yz->ParError(1);
}

GetClusters()

std::vector< Cluster > trigger::IdealClusterBuilder::GetClusters ( )

inline

Definition at line 172 of file IdealClusterBuilder.h.

{ return all_clusters; }

SetClusterGeo()

void trigger::IdealClusterBuilder::SetClusterGeo ( ClusterGeometry *  _g )

inline

Definition at line 173 of file IdealClusterBuilder.h.

{ g = _g; }

Member Data Documentation

all_clusters

std::vector<Cluster> trigger::IdealClusterBuilder::all_clusters {}

Definition at line 146 of file IdealClusterBuilder.h.

{};

all_hits

std::vector<Hit> trigger::IdealClusterBuilder::all_hits {}

Definition at line 145 of file IdealClusterBuilder.h.

{};

debug

bool trigger::IdealClusterBuilder::debug = false

Definition at line 164 of file IdealClusterBuilder.h.

DEPTH_GOOD

int trigger::IdealClusterBuilder::DEPTH_GOOD = 5

Definition at line 160 of file IdealClusterBuilder.h.

g

ClusterGeometry* trigger::IdealClusterBuilder::g

Definition at line 147 of file IdealClusterBuilder.h.

LAYER_MAX

const int trigger::IdealClusterBuilder::LAYER_MAX = 35

Definition at line 155 of file IdealClusterBuilder.h.

LAYER_SEEDMAX

const int trigger::IdealClusterBuilder::LAYER_SEEDMAX = 15

Definition at line 158 of file IdealClusterBuilder.h.

LAYER_SEEDMIN

const int trigger::IdealClusterBuilder::LAYER_SEEDMIN = 3

Definition at line 157 of file IdealClusterBuilder.h.

LAYER_SHOWERMAX

const int trigger::IdealClusterBuilder::LAYER_SHOWERMAX = 7

Definition at line 156 of file IdealClusterBuilder.h.

MIN_TP_ENERGY

const float trigger::IdealClusterBuilder::MIN_TP_ENERGY = 0.5

Definition at line 159 of file IdealClusterBuilder.h.

n_neighbors

int trigger::IdealClusterBuilder::n_neighbors = 1

Definition at line 151 of file IdealClusterBuilder.h.

neighb_thresh

float trigger::IdealClusterBuilder::neighb_thresh = 0

Definition at line 150 of file IdealClusterBuilder.h.

seed_thresh

float trigger::IdealClusterBuilder::seed_thresh = 0

Definition at line 149 of file IdealClusterBuilder.h.

split_energy

bool trigger::IdealClusterBuilder::split_energy = true

Definition at line 152 of file IdealClusterBuilder.h.

use_towers

bool trigger::IdealClusterBuilder::use_towers = false

Definition at line 154 of file IdealClusterBuilder.h.

---

The documentation for this class was generated from the following files:

• Trigger/Algo/include/Trigger/IdealClusterBuilder.h
• Trigger/Algo/src/Trigger/IdealClusterBuilder.cxx