SimChannel.cxx

Go to the documentation of this file.

///
/// \file  Simulation/SimChannel.cxx
///
///
/// \author  seligman@nevis.columbia.edu
///
////////////////////////////////////////////////////////////////////////

#include <limits> // std::numeric_limits
#include <utility>
#include <stdexcept>
#include <algorithm> // std::lower_bound(), std::max()
#include <map>

#include "lardataobj/Simulation/SimChannel.h"
#include "lardataobj/Simulation/sim.h"
#include "larcoreobj/SimpleTypesAndConstants/PhysicalConstants.h"

#include "messagefacility/MessageLogger/MessageLogger.h"

namespace sim{

  //-------------------------------------------------
  IDE::IDE()
    : trackID     (util::kBogusI)
    , numElectrons(util::kBogusD)
    , energy      (util::kBogusD)
    , x           (util::kBogusD)
    , y           (util::kBogusD)
    , z           (util::kBogusD)
  {}

  //-------------------------------------------------
  IDE::IDE(sim::IDE const& ide, int offset)
    : IDE(ide)
  {

    trackID = trackID>=0? trackID+offset : trackID-offset;

  }

  // Default constructor
  //-------------------------------------------------
  SimChannel::SimChannel()
    : fChannel(raw::InvalidChannelID)
  {}

  //-------------------------------------------------
  SimChannel::SimChannel(raw::ChannelID_t channel)
    : fChannel(channel)
  {}

  //-------------------------------------------------
  void SimChannel::AddIonizationElectrons(TrackID_t     trackID,
                                          TDC_t         tdc,
                                          double        numberElectrons,
                                          double const* xyz,
                                          double        energy)
  {
    // look at the collection to see if the current TDC already
    // exists, if not, add it, if so, just add a new track id to the
    // vector, or update the information if track is already present

    // no electrons? no energy? no good!
    if ((numberElectrons < std::numeric_limits<double>::epsilon())
      || (energy <= std::numeric_limits<double>::epsilon()))
    {
      // will throw
      MF_LOG_ERROR("SimChannel")
      << "AddIonizationElectrons() trying to add to TDC #"
      << tdc
      << " "
      << numberElectrons
      << " electrons with "
      << energy
      << " MeV of energy from track ID="
      << trackID;
      return;
    } // if no energy or no electrons

    auto itr = findClosestTDCIDE(tdc);

    // check if this tdc value is in the vector, it is possible that
    // the lower bound is different from the given tdc, in which case
    // we need to add something for that tdc
    if(itr        == fTDCIDEs.end() ||
       itr->first != tdc){
      std::vector<sim::IDE> idelist;
      idelist.emplace_back(trackID,
                           numberElectrons,
                           energy,
                           xyz[0],
                           xyz[1],
                           xyz[2]
                           );
      fTDCIDEs.emplace(itr, tdc, std::move(idelist) );
    }
    else { // we have that TDC already; itr points to it

      // loop over the IDE vector for this tdc and add the electrons
      // to the entry with the same track id
      for(auto& ide : itr->second){

        if (ide.trackID != trackID ) continue;

        // make a weighted average for the location information
        double weight    = ide.numElectrons + numberElectrons;
        ide.x            = (ide.x * ide.numElectrons + xyz[0]*numberElectrons)/weight;
        ide.y            = (ide.y * ide.numElectrons + xyz[1]*numberElectrons)/weight;
        ide.z            = (ide.z * ide.numElectrons + xyz[2]*numberElectrons)/weight;
        ide.numElectrons = weight;
        ide.energy       = ide.energy + energy;

        // found the track id we wanted, so return;
        return;
      } // for

      // if we never found the track id, then this is the first instance of
      // the track id for this tdc, so add ide to the vector
      itr->second.emplace_back(trackID,
                               numberElectrons,
                               energy,
                               xyz[0],
                               xyz[1],
                               xyz[2]
                               );

    } // if new TDC ... else

  } // SimChannel::AddIonizationElectrons()


  //-------------------------------------------------
  double SimChannel::Charge(TDC_t tdc) const
  {
    double charge = 0.;

    auto itr = findClosestTDCIDE(tdc);

      // check to see if this tdc value is in the map
    if(itr        != fTDCIDEs.end() &&
       itr->first == tdc){

        // loop over the list for this tdc value and add up
        // the total number of electrons
      for(auto ide : itr->second){
        charge += ide.numElectrons;
      } // end loop over sim::IDE for this tdc

    } // end if this tdc is represented in the map

    return charge;
  }

    //-------------------------------------------------
  double SimChannel::Energy(TDC_t tdc) const
  {
    double energy = 0.;

    auto itr = findClosestTDCIDE(tdc);

      // check to see if this tdc value is in the map
    if(itr        != fTDCIDEs.end() &&
       itr->first == tdc){

        // loop over the list for this tdc value and add up
        // the total number of electrons
      for(auto ide : itr->second ){
        energy += ide.energy;
      } // end loop over sim::IDE for this tdc

    } // end if this tdc is represented in the map

    return energy;
  }


  //-----------------------------------------------------------------------
  // the start and end tdc values are assumed to be inclusive
  std::vector<sim::IDE> SimChannel::TrackIDsAndEnergies(TDC_t startTDC,
                                                        TDC_t endTDC) const
  {
    // make a map of track ID values to sim::IDE objects

    if(startTDC > endTDC ){
      mf::LogWarning("SimChannel") << "requested tdc range is bogus: "
                                   << startTDC << " " << endTDC
                                   << " return empty vector";
      return {}; // returns an empty vector
    }

    std::map<TrackID_t, sim::IDE> idToIDE;

      //find the lower bound for this tdc and then iterate from there
    auto itr = findClosestTDCIDE(startTDC);

    while(itr != fTDCIDEs.end()){

      // check the tdc value for the iterator, break the loop if we
      // are outside the range
      if(itr->first > endTDC) break;

      // grab the vector of IDEs for this tdc
      auto const& idelist = itr->second;
      // now loop over them and add their content to the map
      for(auto const& ide : idelist){
        auto itTrkIDE = idToIDE.find(ide.trackID);
        if( itTrkIDE != idToIDE.end() ){
          // the IDE we are going to update:
          sim::IDE& trackIDE = itTrkIDE->second;

          double const nel1   = trackIDE.numElectrons;
          double const nel2   = ide.numElectrons;
          double const en1    = trackIDE.energy;
          double const en2    = ide.energy;
          double const energy = en1  + en2;
          double const weight = nel1 + nel2;

            // make a weighted average for the location information
          trackIDE.x            = (ide.x*nel2 + trackIDE.x*nel1)/weight;
          trackIDE.y            = (ide.y*nel2 + trackIDE.y*nel1)/weight;
          trackIDE.z            = (ide.z*nel2 + trackIDE.z*nel1)/weight;
          trackIDE.numElectrons = weight;
          trackIDE.energy = energy;
        } // end if the track id for this one is found
        else{
          idToIDE[ide.trackID] = sim::IDE(ide);
        }
      } // end loop over vector

      ++itr;
    } // end loop over tdc values

      // now fill the vector with the ides from the map
    std::vector<sim::IDE> ides;
    ides.reserve(idToIDE.size());
    for(auto const& itr : idToIDE){
      ides.push_back(itr.second);
    }

    return ides;
  }

  //-----------------------------------------------------------------------
  // the start and end tdc values are assumed to be inclusive
  std::vector<sim::TrackIDE>  SimChannel::TrackIDEs(TDC_t startTDC,
                                                    TDC_t endTDC) const
  {

    std::vector<sim::TrackIDE> trackIDEs;

    if(startTDC > endTDC ){
      mf::LogWarning("SimChannel::TrackIDEs") << "requested tdc range is bogus: "
                                              << startTDC << " " << endTDC
                                              << " return empty vector";
      return trackIDEs;
    }

    double totalE = 0.;
    std::vector<sim::IDE> const ides = TrackIDsAndEnergies(startTDC, endTDC);
    for (auto const& ide : ides)
      totalE += ide.energy;

    // protect against a divide by zero below
    if(totalE < 1.e-5) totalE = 1.;

    // loop over the entries in the map and fill the input vectors
    for (auto const& ide : ides){
      if(ide.trackID == sim::NoParticleId) continue;
      trackIDEs.emplace_back(ide.trackID, ide.energy/totalE, ide.energy, ide.numElectrons);
    }

    return trackIDEs;
  }


  //-----------------------------------------------------------------------
  // Merge the collection of IDEs from one sim channel to another.
  // Requires an agreed upon offset for G4 trackID
  std::pair<SimChannel::TrackID_t,SimChannel::TrackID_t>
  SimChannel::MergeSimChannel(SimChannel const& channel,
                              int               offset)
  {
    if( this->Channel() != channel.Channel() )
      throw std::runtime_error("ERROR SimChannel Merge: Trying to merge different channels!");

    std::pair<TrackID_t,TrackID_t> range_trackID(std::numeric_limits<int>::max(),
                                                 std::numeric_limits<int>::min());

    for(auto const& itr : channel.TDCIDEMap()){

      auto tdc  = itr.first;
      auto const& ides = itr.second;

      // find the entry from this SimChannel corresponding to the tdc from the other
      auto itrthis = findClosestTDCIDE(tdc);

      // pick which IDE list we have to fill: new one or existing one
      std::vector<sim::IDE>* curIDEVec;
      if(itrthis        == fTDCIDEs.end() ||
         itrthis->first != tdc){
        fTDCIDEs.emplace_back(tdc, std::vector<sim::IDE>());
        curIDEVec = &(fTDCIDEs.back().second);
      }
      else
        curIDEVec = &(itrthis->second);

      for(auto const& ide : ides){
        curIDEVec->emplace_back(ide, offset);
        auto tid = std::abs(ide.trackID)+offset;
        if( tid < range_trackID.first  )
          range_trackID.first = tid;
        if( tid > range_trackID.second )
          range_trackID.second = tid;
      }//end loop over IDEs

    }//end loop over TDCIDEMap


    return range_trackID;

  }

  //-------------------------------------------------
  struct SimChannel::CompareByTDC {

    bool operator()
      (TDCIDE const& a, TDCIDE const& b) const
      { return a.first < b.first; }

    bool operator()
      (StoredTDC_t a_tdc, TDCIDE const& b) const
      { return a_tdc < b.first; }

    bool operator()
      (TDCIDE const& a, StoredTDC_t b_tdc) const
      { return a.first < b_tdc; }

  }; // struct CompareByTDC


  SimChannel::TDCIDEs_t::iterator SimChannel::findClosestTDCIDE(StoredTDC_t tdc)
  {
    return std::lower_bound
      (fTDCIDEs.begin(), fTDCIDEs.end(), tdc, CompareByTDC());
  }

  SimChannel::TDCIDEs_t::const_iterator SimChannel::findClosestTDCIDE
    (StoredTDC_t tdc) const
  {
    return std::lower_bound
      (fTDCIDEs.begin(), fTDCIDEs.end(), tdc, CompareByTDC());
  }


  //-------------------------------------------------


}