TrigScintClusterProducer.cxx

 
#include "TrigScint/TrigScintClusterProducer.h"
 
#include <iterator>  // std::next
#include <map>
 
namespace trigscint {
 
void TrigScintClusterProducer::configure(framework::config::Parameters &ps) {
  seed_ = ps.get<double>("seed_threshold");
  min_thr_ = ps.get<double>("clustering_threshold");
  max_width_ = ps.get<int>("max_cluster_width");
  ampl_weighting_ = ps.get<bool>("ampl_weighting");
  input_collection_ = ps.get<std::string>("input_collection");
  pass_name_ = ps.get<std::string>("input_pass_name");
  output_collection_ = ps.get<std::string>("output_collection");
  verbose_ = ps.get<int>("verbosity");
  vert_bar_start_idx_ = ps.get<int>("vertical_bar_start_index");
  time_tolerance_ = ps.get<double>("time_tolerance");
  pad_time_ = ps.get<double>("pad_time");
  if (verbose_) {
    ldmx_log(info) << "In TrigScintClusterProducer: configure done!";
    ldmx_log(info) << "Got parameters: \nSeed threshold:   " << seed_
                   << "\nClustering threshold: " << min_thr_
                   << "\nMax cluster width: " << max_width_
                   << "\nAmplitude weighting: " << ampl_weighting_
                   << "\nExpected pad hit time: " << pad_time_
                   << "\nMax hit time delay: " << time_tolerance_
                   << "\nVertical bar start index:     " << vert_bar_start_idx_
                   << "\nInput collection:     " << input_collection_
                   << "\nInput pass name:     " << pass_name_
                   << "\nOutput collection:    " << output_collection_
                   << "\nVerbosity: " << verbose_;
  }
 
  return;
}
 
void TrigScintClusterProducer::produce(framework::Event &event) {
  // parameters.
  // a cluster seeding threshold
  // a clustering threshold -- a lower boundary for being added at all (zero
  // suppression)
  // a maximum cluster width
  //
  // steps.
  // 1. get an input collection of digi hits. at most one entry per channel.
  // 2. access them by channel number
  // 3. clustering:
  //       a. add first hit > seedThr to cluster(ling) . store content as
  //       localMax. track size of cluster (1) b. if not in beginning of array:
  //       add cell before while content <  add next hit first hit > seedThr to
  //       cluster(ling) b. while content <  add next hit first hit > seedThr to
  //       cluster(ling)
 
  /*
 
    The trigger pad geometry considered for this clustering algorithm is:
 
 
    |   |   |   |   |   |
    | 0 | 2 | 4 | 6 | 8 |  ... | 48|
 
      | 1 | 3 | 5 | 7 | 9 |  ... | 49|
      |   |   |   |   |   |
 
    with hits in channels after digi looking something like this
 
 
    ampl:    _                   _                   _
            | |_                | |                 | |
        ----| | |---------------| |-----------------| |------- cluster seed
    threshold | | |_              | |_ _              | |_
           _| | | |            _| | | |            _| | |  _
          | | | | |     vs    | | | | |    vs     | | | | | |
 
 
              |                   |                   |
        split | seeds        keep | disregard   keep! | just move on
              | next cl.          | (later seed       | (no explicit splitting)
                                  might pick
                                  it up)
 
 
    The idea being that while there could be good reasons for an electron to
    touch three pads in a row, there is no good reason for it to cross four.
    This is noise, or, the start of an adjacent cluster. In any case, 4 is not a
    healthy cluster. Proximity to a seed governs which below-seed channels to
    include. By always starting in one end but going back (at most two
    channels), this algo guarantees symmetric treatment on both sides of the
    seed.
 
 
    //Procedure: keep going until there is a seed. walk back at most 2 steps
    // add all the hits. clusters of up to 3 is fine.
    // if the cluster is > 3, then we need to do something.
    // if it's == 4, we'd want to split in the middle if there are two potential
    seeds. retain only the first half, cases are (seed = s, n - no/noise)
    // nsns , nssn, snsn, ssnn.  nnss won't happen (max 1 step back from s,
    unless there is nothing in front)
    // all these are ok to split like this bcs even if ssnn--> ss and some small
    nPE is lost, that's probably pretty negligible wrt the centroid position,
    with two seeds in one cluster
 
    // if it's > 4, cases are
    // nsnsn, nsnss, nssnn, nssns, snsnn, snsns, ssnnn, ssnns.
    // these are also all ok to just truncate after 2. and then the same check
    outlined above will happen to the next chunk.
 
    // so in short we can
    // 1. seed --> addHit
    // 2. walk back once --> addHit
    // 3. check next: if seed+1 exists && seed +2 exists,
    // 3a. if seed-1 is in already, stop here.
    // 3b. else if seed+3 exists, stop here.
    // 3c. else addHit(seed+1), addHit(seed+2)
    // 4. if seed+1 and !seed+2 --> addHit(seed+1)
    // 5. at this point, if clusterSize is 2 hits and seed+1 didn't exist, we
    can afford to walk back one more step and add whatever junk was there (we
    know it's not a seed)
 
 
    */
 
  if (verbose_) {
    ldmx_log(debug)
        << "TrigScintClusterProducer: produce() starts! Event number: "
        << event.getEventHeader().getEventNumber();
  }
 
  // looper over digi hits and aggregate energy depositions for each detID
 
  const auto digis{
      event.getCollection<ldmx::TrigScintHit>(input_collection_, pass_name_)};
 
  if (verbose_) {
    ldmx_log(debug) << "Got digi collection " << input_collection_ << "_"
                    << pass_name_ << " with " << digis.size() << " entries ";
  }
 
  // TODO remove this once verified that the noise overlap bug is gone
  bool do_duplicate = true;
 
  // 1. store all the channel digi content in channel order
  auto i_digi{0};
  for (const auto &digi : digis) {
    // these are unordered hits, and this collection is zero-suppressed
    // map the index of the digi to the channel index
 
    if (digi.getPE() >
        min_thr_) {  // cut on a min threshold (for a non-seeding hit to be
                     // added to seeded clusters) already here
 
      int id = digi.getBarID();
 
      // first check if there is a (pure) noise hit at this channel,  and
      // replace it in that case. this is a protection against a problem that
      // shouldn't be there in the first place.
      if (do_duplicate &&
          hit_channel_map_.find((id)) != hit_channel_map_.end()) {
        int idx = id;
        std::map<int, int>::iterator itr = hit_channel_map_.find(idx);
        double old_val = digis.at(itr->second).getPE();
        if (verbose_) {
          ldmx_log(debug) << "Got duplicate digis for channel " << idx
                          << ", with already inserted value " << old_val
                          << " and new " << digi.getPE();
        }
        if (digi.getPE() > old_val) {
          hit_channel_map_.erase(itr->first);
          if (verbose_) {
            ldmx_log(debug)
                << "Skipped duplicate digi with smaller value for channel "
                << idx;
          }
        }
      }
 
      // don't add in late hits
      if (digi.getTime() > pad_time_ + time_tolerance_) {
        i_digi++;
        continue;
      }
 
      hit_channel_map_.insert(std::pair<int, int>(id, i_digi));
      // the channel number is the key, the digi list index is the value
 
      if (verbose_) {
        ldmx_log(debug) << "Mapping digi hit nb " << i_digi
                        << " with energy = " << digi.getEnergy()
                        << " MeV, nPE = " << digi.getPE() << " > " << min_thr_
                        << " to key/channel " << id;
      }
    }
    i_digi++;
  }
 
  // 2. now step through all the channels in the map and cluster the hits
 
  std::map<int, int>::iterator itr;
 
  // Create the container to hold the digitized trigger scintillator hits.
  std::vector<ldmx::TrigScintCluster> trig_scint_clusters;
 
  // loop over channels
  for (itr = hit_channel_map_.begin(); itr != hit_channel_map_.end(); ++itr) {
    // this hit may have disappeared
    if (hit_channel_map_.find(itr->first) == hit_channel_map_.end()) {
      if (verbose_ > 1) {
        ldmx_log(debug) << "Attempting to use removed hit at channel "
                        << itr->first << "; skipping.";
      }
      continue;
    }
 
    // i don't like this but for now, erasing elements in the map leads, as it
    // turns out, to edge cases where i miss out on hits or run into
    // non-existing indices. so while what i do below means that i don't need to
    // erase hits, i'd rather find a way to do that and skip this book keeping:
    bool has_used = false;
    for (const auto &index : v_used_indices_) {
      if (index == itr->first) {
        if (verbose_ > 1) {
          ldmx_log(warn) << "Attempting to re-use hit at channel " << itr->first
                         << "; skipping.";
        }
        has_used = true;
      }
    }
    if (has_used) continue;
    if (verbose_ > 1) {
      ldmx_log(debug) << "\t At hit with channel nb " << itr->first << ".";
    }
 
    if (hit_channel_map_.size() ==
        0)  // we removed them all..? shouldn't ever happen
    {
      if (verbose_)
        ldmx_log(warn) << "Time flies, and all clusters have already been "
                          "removed! Unclear how we even got here; interfering "
                          "here to get out of the loop. ";
      break;
    }
 
    ldmx::TrigScintHit digi = (ldmx::TrigScintHit)digis.at(itr->second);
 
    // skip all until hit a seed
    if (digi.getPE() >= seed_) {
      if (verbose_ > 1) {
        ldmx_log(debug) << "Seeding cluster with channel " << itr->first
                        << "; content " << digi.getPE();
      }
 
      // 1.  add seeding hit to cluster
 
      addHit(itr->first, digi);
 
      if (verbose_ > 1) {
        ldmx_log(debug) << "\t itr is pointing at hit with channel nb "
                        << itr->first << ".";
      }
 
      // ----- first look back one step
 
      // we have added the hit from the neighbouring channel to the list only if
      // it's above clustering threshold so no check needed now
      std::map<int, int>::iterator itr_back =
          hit_channel_map_.find(itr->first - 1);
 
      bool has_backed = false;
 
      if (itr_back !=
          hit_channel_map_
              .end()) {  // there is an entry for the previous
                         // channel, so it had content above threshold
        // but it wasn't enough to seed a cluster. so, unambiguous that it
        // should be added here because it's its only chance to get in.
 
        // need to check again for backwards hits
        has_used = false;
        for (const auto &index : v_used_indices_) {
          if (index == itr_back->first) {
            if (verbose_ > 1) {
              ldmx_log(warn) << "Attempting to re-use hit at channel "
                             << itr_back->first << "; skipping.";
            }
            has_used = true;
          }
        }
        if (!has_used) {
          digi = (ldmx::TrigScintHit)digis.at(itr_back->second);
 
          // 2. add seed-1 to cluster
          addHit(itr_back->first, digi);
          has_backed = true;
 
          if (verbose_ > 1) {
            ldmx_log(debug) << "Added -1 channel " << itr_back->first
                            << " to cluster; content " << digi.getPE();
            ldmx_log(debug) << "\t itr is pointing at hit with channel nb "
                            << itr->first << ".";
          }
 
        }  // if seed-1 wasn't used already
      }  // there exists a lower, unused neighbour
 
      // 3. check next: if seed+1 exists && seed +2 exists,
      //    3a. if seed-1 is in already, this is a case for a split, at seed. go
      //    directly to check on seed-2, don't add more here. 3b. else. addHit
      //    (seed+1) 3c. if seed+3 exists, this is a split, at seed+1. don't add
      //    more here. 3d. else addHit(seed+2)
      // 4. if seed+1 and !seed+2 --> go to addHit(seed+1)
 
      // --- now, step 3, 4: look ahead 1 step from seed
 
      if (v_added_indices_.size() < max_width_) {
        // (in principle these don't need to be different iterators, but it
        // makes the logic easier to follow)
        std::map<int, int>::iterator itr_neighb =
            hit_channel_map_.find(itr->first + 1);
        if (itr_neighb !=
            hit_channel_map_
                .end()) {  // there is an entry for the next channel,
                           // so it had content above threshold
          // seed+1 exists
          // check if there is sth in position seed+2
          if (hit_channel_map_.find(itr_neighb->first + 1) !=
              hit_channel_map_.end()) {  // a hit with that key exists, so
                                         // seed+1 and seed+2 exist
            if (!has_backed) {  // there is no seed-1 in the cluster. room for
                                // at least seed+1, and for seed+2 only if there
                                // is no seed+3
              // 3b
              digi = (ldmx::TrigScintHit)digis.at(itr_neighb->second);
              addHit(itr_neighb->first, digi);
 
              if (verbose_ > 1) {
                ldmx_log(debug)
                    << "No -1 hit. Added +1 channel " << itr_neighb->first
                    << " to cluster; content " << digi.getPE();
                ldmx_log(debug) << "\t itr is pointing at hit with channel nb "
                                << itr->first << ".";
              }
 
              if (v_added_indices_.size() < max_width_) {
                if (hit_channel_map_.find(itr_neighb->first + 2) ==
                    hit_channel_map_
                        .end()) {  // no seed+3. also no seed-1. so add seed+2
                  // 3d.  add seed+2 to the cluster
                  itr_neighb = hit_channel_map_.find(itr->first + 2);
                  digi = (ldmx::TrigScintHit)digis.at(itr_neighb->second);
                  addHit(itr_neighb->first, digi);
                  if (verbose_ > 1) {
                    ldmx_log(debug)
                        << "No +3 hit. Added +2 channel " << itr_neighb->first
                        << " to cluster; content " << digi.getPE();
                    ldmx_log(debug)
                        << "\t itr is pointing at hit with channel nb "
                        << itr->first << ".";
                  }
                }
 
              }  // if no seed+3 --> added seed+2
            }  // if seed-1 wasn't added
          }  // if seed+2 exists. then already added seed+1.
          else {  // so: if not, then we need to add seed+1 here. (step 4)
            digi = (ldmx::TrigScintHit)digis.at(
                itr_neighb->second);  // itrNeighb hasn't moved since there was
                                      // no seed+2
            addHit(itr_neighb->first, digi);
 
            if (verbose_ > 1) {
              ldmx_log(debug)
                  << "Added +1 channel " << itr_neighb->first
                  << " as last channel to cluster; content " << digi.getPE();
              ldmx_log(debug) << "\t itr is pointing at hit with channel nb "
                              << itr->first << ".";
            }
          }
        }  // if seed+1 exists
        // 5. at this point, if clusterSize is 2 hits and seed+1 didn't exist,
        // we can afford to walk back one more step and add whatever junk was
        // there (we know it's not a seed)
        else if (has_backed &&
                 hit_channel_map_.find(itr_back->first - 1) !=
                     hit_channel_map_
                         .end()) {  // seed-1 has been added, but not seed+1,
                                    // and there is a hit in seed-2
          itr_back = hit_channel_map_.find(itr->first - 2);
          digi = (ldmx::TrigScintHit)digis.at(itr_back->second);
          addHit(itr_back->first, digi);
 
          if (verbose_ > 1) {
            ldmx_log(debug) << "Added -2 channel " << itr_back->first
                            << " to cluster; content " << digi.getPE();
          }
          if (verbose_ > 1) {
            ldmx_log(debug) << "\t itr is pointing at hit with channel nb "
                            << itr->first << ".";
          }
 
        }  // check if add in seed -2
 
      }  // if adding another hit, going forward, was allowed
 
      // done adding hits to cluster. calculate centroid
      centroid_ /=
          sumw_;  // final weighting step: divide by total amplitude sum
 
      centroid_ -= 1;  // shift back to actual channel center
 
      ldmx::TrigScintCluster cluster;
 
      if (verbose_ > 1) {
        ldmx_log(debug) << "Now have " << v_added_indices_.size()
                        << " hits in the cluster ";
      }
      cluster.setSeed(v_added_indices_.at(0));
      cluster.setIDs(v_added_indices_);
      cluster.setNHits(v_added_indices_.size());
      cluster.setCentroid(centroid_);
      float cx;
      float cy = centroid_;
      float cz = -99999;  // set to nonsense for now. could be set to module nb
      // then in horizontal bars --> we don't know X
      if (centroid_ < vert_bar_start_idx_) {
        // set to nonsense in barID space. could translate to x=0 mm
        cx = -1;
      }
 
      else {
        cx = (int)((centroid_ - vert_bar_start_idx_) / 4);  // start at 0
        cy = (int)centroid_ % 4;
      }
      cluster.setCentroidXYZ(cx, cy, cz);
      cluster.setEnergy(val_e_);
      cluster.setPE(val_);
      cluster.setTime(time_ / val_);
      cluster.setBeamEfrac(beam_e_ / val_e_);
 
      trig_scint_clusters.push_back(cluster);
 
      ldmx_log(trace) << cluster;
 
      centroid_ = 0;
      centroid_x_ = -1;
      centroid_y_ = -1;
      val_ = 0;
      val_e_ = 0;
      beam_e_ = 0;
      time_ = 0;
      sumw_ = 0;
      // book keep which channels have already been added to a cluster
      v_added_indices_.resize(0);
 
      if (verbose_ > 1) {
        ldmx_log(debug)
            << "\t Finished processing of seeding hit with channel nb "
            << itr->first << ".";
      }
 
    }  // if content enough to seed a cluster
 
    if (hit_channel_map_.begin() == hit_channel_map_.end()) {
      if (verbose_)
        ldmx_log(warn) << "Time flies, and all clusters have already been "
                          "removed! Interfering here to get out of the loop. ";
      break;
    }
  }  // over channels
 
  if (trig_scint_clusters.size() > 0)
    event.add(output_collection_, trig_scint_clusters);
 
  hit_channel_map_.clear();
  // book keep which channels have already been added to a cluster
  v_used_indices_.resize(0);
 
  return;
}
 
void TrigScintClusterProducer::addHit(uint idx, ldmx::TrigScintHit hit) {
  float ampl = hit.getPE();
  float w = 1;
  if (ampl_weighting_) {  // if choosing to PE-weight centroid positions
    w = ampl;
  }
 
  float energy = hit.getEnergy();
  val_e_ += energy;
 
  val_ += ampl;
  centroid_ += (idx + 1) * w;  // need non-zero weight of channel 0. shifting
                               // centroid back by 1 in the end
                               // this number gets divided by val at the end
 
  sumw_ += w;
 
  v_added_indices_.push_back(idx);
 
  beam_e_ += hit.getBeamEfrac() * energy;
  if (hit.getTime() > -990.) {
    time_ += hit.getTime() * ampl;
  }
 
  v_used_indices_.push_back(idx);
  /*    // not working properly, but i'd prefer this type of solution
        hit_channel_map_.erase( idx ) ;
        if (verbose_ > 1 ) {
        ldmx_log(debug) << "Removed used hit " << idx << " from list";
        }
        if ( hit_channel_map_.find( idx) != hit_channel_map_.end() )
        std::cerr << "----- WARNING! Hit still present in map after removal!! ";
  */
  if (verbose_ > 1) {
    ldmx_log(debug) << "   In addHit, adding hit at " << idx
                    << " with amplitude " << ampl
                    << ", updating cluster to current centroid "
                    << centroid_ / val_ - 1 << " and energy " << val_
                    << ". index vector now ends with "
                    << v_added_indices_.back();
  }
 
  return;
}
 
void TrigScintClusterProducer::onProcessStart() {
  ldmx_log(debug) << "Process starts!";
 
  return;
}
 
void TrigScintClusterProducer::onProcessEnd() {
  ldmx_log(debug) << "Process ends!";
 
  return;
}
 
}  // namespace trigscint
 
DECLARE_PRODUCER(trigscint::TrigScintClusterProducer);