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
const int	DEPTH_GOOD = 5
bool	debug_ = false

Detailed Description

Definition at line 140 of file IdealClusterBuilder.h.

Member Function Documentation

addHit()

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

inline

Definition at line 164 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 193 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_) {
    layer_hits[hit.layer_].push_back(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 49 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 is_local_max = true;
    for (auto n : g_->neighbors_[hit.id_]) {
      if (hits_by_id.count(n) && hits_by_id[n].e_ > hit.e_)
        is_local_max = false;
    }
    if (is_local_max && (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_clus2hit_i_ds;
    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_clus2hit_i_ds[iclus] = neighbors;
    }
 
    // check how many clusters to which each hit is assoc
    std::map<int, std::vector<int> > assoc_hit_i_d2clusters;
    for (auto clus2hit_id : assoc_clus2hit_i_ds) {
      auto iclus = clus2hit_id.first;
      auto& hit_i_ds = clus2hit_id.second;
      for (const auto& hit_id : hit_i_ds) {
        assoc_hit_i_d2clusters[hit_id].push_back(iclus);
      }
    }
 
    // add associated hits to clusters
    //   (w/ optional e-splitting)
    for (auto hit_i_d2clusters : assoc_hit_i_d2clusters) {
      auto hit_id = hit_i_d2clusters.first;
      auto iclusters = hit_i_d2clusters.second;
      if (iclusters.size() == 1) {
        // simply add cell to the cluster
        auto& hit = hits_by_id[hit_id];
        auto iclus = iclusters[0];
        hit.used_ = true;
        clusters[iclus].hits_.push_back(hit);
        clusters[iclus].e_ += hit.e_;
      } else {
        auto& hit = hits_by_id[hit_id];
        hit.used_ = true;
        float esum = 0;
        for (auto iclus : iclusters) {
          esum += clusters[iclus].e_;
        }
        for (auto iclus : iclusters) {
          Hit new_hit = hit;
          if (split_energy_) new_hit.e_ = hit.e_ * clusters[iclus].e_ / esum;
          clusters[iclus].hits_.push_back(new_hit);
          clusters[iclus].e_ += new_hit.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 211 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.is_2d_ = 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_.clear();
          cluster3d.clusters2d_.push_back(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 394 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_].n_sub_hit_++;
      } else {
        towers[hit.id_] = hit;
        towers[hit.id_].layer_ = 0;
        towers[hit.id_].z_ = 0;
        towers[hit.id_].n_sub_hit_ = 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 432 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 170 of file IdealClusterBuilder.h.

{ return all_clusters_; }

setClusterGeo()

void trigger::IdealClusterBuilder::setClusterGeo ( ClusterGeometry * g )

inline

Definition at line 171 of file IdealClusterBuilder.h.

{ g_ = g; }

Member Data Documentation

all_clusters_

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

Definition at line 144 of file IdealClusterBuilder.h.

{};

all_hits_

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

Definition at line 143 of file IdealClusterBuilder.h.

{};

debug_

bool trigger::IdealClusterBuilder::debug_ = false

Definition at line 162 of file IdealClusterBuilder.h.

DEPTH_GOOD

const int trigger::IdealClusterBuilder::DEPTH_GOOD = 5

Definition at line 158 of file IdealClusterBuilder.h.

g_

ClusterGeometry* trigger::IdealClusterBuilder::g_

Definition at line 145 of file IdealClusterBuilder.h.

LAYER_MAX

const int trigger::IdealClusterBuilder::LAYER_MAX = 35

Definition at line 153 of file IdealClusterBuilder.h.

LAYER_SEEDMAX

const int trigger::IdealClusterBuilder::LAYER_SEEDMAX = 15

Definition at line 156 of file IdealClusterBuilder.h.

LAYER_SEEDMIN

const int trigger::IdealClusterBuilder::LAYER_SEEDMIN = 3

Definition at line 155 of file IdealClusterBuilder.h.

LAYER_SHOWERMAX

const int trigger::IdealClusterBuilder::LAYER_SHOWERMAX = 7

Definition at line 154 of file IdealClusterBuilder.h.

MIN_TP_ENERGY

const float trigger::IdealClusterBuilder::MIN_TP_ENERGY = 0.5

Definition at line 157 of file IdealClusterBuilder.h.

n_neighbors_

int trigger::IdealClusterBuilder::n_neighbors_ = 1

Definition at line 149 of file IdealClusterBuilder.h.

neighb_thresh_

float trigger::IdealClusterBuilder::neighb_thresh_ = 0

Definition at line 148 of file IdealClusterBuilder.h.

seed_thresh_

float trigger::IdealClusterBuilder::seed_thresh_ = 0

Definition at line 147 of file IdealClusterBuilder.h.

split_energy_

bool trigger::IdealClusterBuilder::split_energy_ = true

Definition at line 150 of file IdealClusterBuilder.h.

use_towers_

bool trigger::IdealClusterBuilder::use_towers_ = false

Definition at line 152 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