PHG4MicromegasDigitizer.cc

Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHG4MicromegasDigitizer.cc

#include "PHG4MicromegasDigitizer.h"

// Move to new storage containers
#include <trackbase/TrkrDefs.h>
#include <trackbase/TrkrHit.h>
#include <trackbase/TrkrHitSet.h>
#include <trackbase/TrkrHitSetContainer.h>
#include <trackbase/TrkrHitTruthAssoc.h>

#include <phparameter/PHParameterInterface.h>  // for PHParameterInterface

#include <fun4all/Fun4AllReturnCodes.h>
#include <fun4all/SubsysReco.h>

#include <phool/PHRandomSeed.h>
#include <phool/getClass.h>

#include <gsl/gsl_randist.h>
#include <gsl/gsl_rng.h>  // for gsl_rng_alloc, gsl_rng...

#include <cassert>
#include <iostream>  // for operator<<, basic_ostream
#include <set>
#include <utility>  // for pair

namespace
{

  // local version of std::clamp, which is only available for c++17
  template <class T>
  constexpr const T& clamp(const T& v, const T& lo, const T& hi)
  {
    return (v < lo) ? lo : (hi < v) ? hi
                                    : v;
  }

}  // namespace
//____________________________________________________________________________
PHG4MicromegasDigitizer::PHG4MicromegasDigitizer(const std::string& name)
  : SubsysReco(name)
  , PHParameterInterface(name)
{
  // initialize rng
  const unsigned int seed = PHRandomSeed();
  m_rng.reset(gsl_rng_alloc(gsl_rng_mt19937));
  gsl_rng_set(m_rng.get(), seed);

  InitializeParameters();
}

//____________________________________________________________________________
int PHG4MicromegasDigitizer::InitRun(PHCompositeNode* /*topNode*/)
{
  UpdateParametersWithMacro();

  // load parameters
  m_adc_threshold = get_double_param("micromegas_adc_threshold");
  m_enc = get_double_param("micromegas_enc");
  m_pedestal = get_double_param("micromegas_pedestal");
  m_volts_per_charge = get_double_param("micromegas_volts_per_charge");

  // printout
  std::cout
      << "PHG4MicromegasDigitizer::InitRun\n"
      << " m_adc_threshold: " << m_adc_threshold << " electrons\n"
      << " m_enc: " << m_enc << " electrons\n"
      << " m_pedestal: " << m_pedestal << " electrons\n"
      << " m_volts_per_charge: " << m_volts_per_charge << " mV/fC\n"
      << std::endl;

  // threshold is effectively applied on top of pedestal
  m_adc_threshold += m_pedestal;

  /*
   * Factor that convertes charge in a voltage in each z bin
   * the scale up factor of 2.4 is meant to account for shaping time (80ns)
   * it only applies to the signal
   * see: simulations/g4simulations/g4tpc/PHG4TpcDigitizer::InitRun
   */
  m_volt_per_electron_signal = m_volts_per_charge * 1.602e-4 * 2.4;
  m_volt_per_electron_noise = m_volts_per_charge * 1.602e-4;

  return Fun4AllReturnCodes::EVENT_OK;
}

//____________________________________________________________________________
int PHG4MicromegasDigitizer::process_event(PHCompositeNode* topNode)
{
  // Get Nodes
  // Get the TrkrHitSetContainer node
  auto trkrhitsetcontainer = findNode::getClass<TrkrHitSetContainer>(topNode, "TRKR_HITSET");
  assert(trkrhitsetcontainer);

  // Get the TrkrHitTruthAssoc node
  auto hittruthassoc = findNode::getClass<TrkrHitTruthAssoc>(topNode, "TRKR_HITTRUTHASSOC");

  // get all micromegas hitsets
  const auto hitset_range = trkrhitsetcontainer->getHitSets(TrkrDefs::TrkrId::micromegasId);
  for (auto hitset_it = hitset_range.first; hitset_it != hitset_range.second; ++hitset_it)
  {
    // get key
    const TrkrDefs::hitsetkey hitsetkey = hitset_it->first;

    // get all of the hits from this hitset
    TrkrHitSet* hitset = hitset_it->second;
    TrkrHitSet::ConstRange hit_range = hitset->getHits();

    // keep track of hits to be removed
    std::set<TrkrDefs::hitkey> removed_keys;

    // loop over hits
    for (auto hit_it = hit_range.first; hit_it != hit_range.second; ++hit_it)
    {
      // store key and hit
      const TrkrDefs::hitkey& key = hit_it->first;
      TrkrHit* hit = hit_it->second;

      // get energy (electrons)
      const double signal = hit->getEnergy();
      const double noise = add_noise();

      // convert to mV
      const double voltage = (m_pedestal + noise) * m_volt_per_electron_noise + signal * m_volt_per_electron_signal;

      // compare to threshold
      if (voltage > m_adc_threshold * m_volt_per_electron_noise)
      {
        // keep hit, update adc
        hit->setAdc(clamp<uint>(voltage * m_adc_per_volt, 0, 1023));
      }
      else
      {
        // mark hit as removable
        removed_keys.insert(key);
      }
    }

    // remove hits
    for (const auto& key : removed_keys)
    {
      hitset->removeHit(key);
      if (hittruthassoc)
      {
        hittruthassoc->removeAssoc(hitsetkey, key);
      }
    }
  }
  return Fun4AllReturnCodes::EVENT_OK;
}

//___________________________________________________________________________
void PHG4MicromegasDigitizer::SetDefaultParameters()
{
  // all values taken from TPC sampa chips (simulations/g4simulations/g4tpc/PHG4TpcDigitizer)
  set_default_double_param("micromegas_enc", 670);
  set_default_double_param("micromegas_adc_threshold", 2680);
  set_default_double_param("micromegas_pedestal", 50000);
  set_default_double_param("micromegas_volts_per_charge", 20);
}

//___________________________________________________________________________
double PHG4MicromegasDigitizer::add_noise() const
{
  return gsl_ran_gaussian(m_rng.get(), m_enc);
}