g4simulation/g4micromegas/PHG4MicromegasDigitizer.cc

0001 /*!
0002  * \file PHG4MicromegasDigitizer.cc
0003  * \author Hugo Pereira Da Costa <hugo.pereira-da-costa@cea.fr>
0004  */
0005
0006 #include "PHG4MicromegasDigitizer.h"
0007
0008 // Move to new storage containers
0009 #include <trackbase/TrkrDefs.h>
0010 #include <trackbase/TrkrHit.h>
0011 #include <trackbase/TrkrHitSet.h>
0012 #include <trackbase/TrkrHitSetContainer.h>
0013 #include <trackbase/TrkrHitTruthAssoc.h>
0014
0015 #include <phparameter/PHParameterInterface.h>  // for PHParameterInterface
0016
0017 #include <fun4all/Fun4AllReturnCodes.h>
0018 #include <fun4all/SubsysReco.h>
0019
0020 #include <phool/PHRandomSeed.h>
0021 #include <phool/getClass.h>
0022
0023 #include <gsl/gsl_randist.h>
0024 #include <gsl/gsl_rng.h>  // for gsl_rng_alloc, gsl_rng...
0025
0026 #include <cassert>
0027 #include <iostream>  // for operator<<, basic_ostream
0028 #include <set>
0029 #include <utility>  // for pair
0030
0031 namespace
0032 {
0033
0034   // local version of std::clamp, which is only available for c++17
0035   template <class T>
0036   constexpr const T& clamp(const T& v, const T& lo, const T& hi)
0037   {
0038     return (v < lo) ? lo : (hi < v) ? hi
0039                                     : v;
0040   }
0041
0042 }  // namespace
0043 //____________________________________________________________________________
0044 PHG4MicromegasDigitizer::PHG4MicromegasDigitizer(const std::string& name)
0045   : SubsysReco(name)
0046   , PHParameterInterface(name)
0047 {
0048   // initialize rng
0049   const unsigned int seed = PHRandomSeed();
0050   m_rng.reset(gsl_rng_alloc(gsl_rng_mt19937));
0051   gsl_rng_set(m_rng.get(), seed);
0052
0053   InitializeParameters();
0054 }
0055
0056 //____________________________________________________________________________
0057 int PHG4MicromegasDigitizer::InitRun(PHCompositeNode* /*topNode*/)
0058 {
0059   UpdateParametersWithMacro();
0060
0061   // load parameters
0062   m_adc_threshold = get_double_param("micromegas_adc_threshold");
0063   m_enc = get_double_param("micromegas_enc");
0064   m_pedestal = get_double_param("micromegas_pedestal");
0065   m_volts_per_charge = get_double_param("micromegas_volts_per_charge");
0066
0067   // printout
0068   std::cout
0069       << "PHG4MicromegasDigitizer::InitRun\n"
0070       << " m_adc_threshold: " << m_adc_threshold << " electrons\n"
0071       << " m_enc: " << m_enc << " electrons\n"
0072       << " m_pedestal: " << m_pedestal << " electrons\n"
0073       << " m_volts_per_charge: " << m_volts_per_charge << " mV/fC\n"
0074       << std::endl;
0075
0076   // threshold is effectively applied on top of pedestal
0077   m_adc_threshold += m_pedestal;
0078
0079   /*
0080    * Factor that convertes charge in a voltage in each z bin
0081    * the scale up factor of 2.4 is meant to account for shaping time (80ns)
0082    * it only applies to the signal
0083    * see: simulations/g4simulations/g4tpc/PHG4TpcDigitizer::InitRun
0084    */
0085   m_volt_per_electron_signal = m_volts_per_charge * 1.602e-4 * 2.4;
0086   m_volt_per_electron_noise = m_volts_per_charge * 1.602e-4;
0087
0088   return Fun4AllReturnCodes::EVENT_OK;
0089 }
0090
0091 //____________________________________________________________________________
0092 int PHG4MicromegasDigitizer::process_event(PHCompositeNode* topNode)
0093 {
0094   // Get Nodes
0095   // Get the TrkrHitSetContainer node
0096   auto trkrhitsetcontainer = findNode::getClass<TrkrHitSetContainer>(topNode, "TRKR_HITSET");
0097   assert(trkrhitsetcontainer);
0098
0099   // Get the TrkrHitTruthAssoc node
0100   auto hittruthassoc = findNode::getClass<TrkrHitTruthAssoc>(topNode, "TRKR_HITTRUTHASSOC");
0101
0102   // get all micromegas hitsets
0103   const auto hitset_range = trkrhitsetcontainer->getHitSets(TrkrDefs::TrkrId::micromegasId);
0104   for (auto hitset_it = hitset_range.first; hitset_it != hitset_range.second; ++hitset_it)
0105   {
0106     // get key
0107     const TrkrDefs::hitsetkey hitsetkey = hitset_it->first;
0108
0109     // get all of the hits from this hitset
0110     TrkrHitSet* hitset = hitset_it->second;
0111     TrkrHitSet::ConstRange hit_range = hitset->getHits();
0112
0113     // keep track of hits to be removed
0114     std::set<TrkrDefs::hitkey> removed_keys;
0115
0116     // loop over hits
0117     for (auto hit_it = hit_range.first; hit_it != hit_range.second; ++hit_it)
0118     {
0119       // store key and hit
0120       const TrkrDefs::hitkey& key = hit_it->first;
0121       TrkrHit* hit = hit_it->second;
0122
0123       // get energy (electrons)
0124       const double signal = hit->getEnergy();
0125       const double noise = add_noise();
0126
0127       // convert to mV
0128       const double voltage = (m_pedestal + noise) * m_volt_per_electron_noise + signal * m_volt_per_electron_signal;
0129
0130       // compare to threshold
0131       if (voltage > m_adc_threshold * m_volt_per_electron_noise)
0132       {
0133         // keep hit, update adc
0134         hit->setAdc(clamp<uint>(voltage * m_adc_per_volt, 0, 1023));
0135       }
0136       else
0137       {
0138         // mark hit as removable
0139         removed_keys.insert(key);
0140       }
0141     }
0142
0143     // remove hits
0144     for (const auto& key : removed_keys)
0145     {
0146       hitset->removeHit(key);
0147       if (hittruthassoc)
0148       {
0149         hittruthassoc->removeAssoc(hitsetkey, key);
0150       }
0151     }
0152   }
0153   return Fun4AllReturnCodes::EVENT_OK;
0154 }
0155
0156 //___________________________________________________________________________
0157 void PHG4MicromegasDigitizer::SetDefaultParameters()
0158 {
0159   // all values taken from TPC sampa chips (simulations/g4simulations/g4tpc/PHG4TpcDigitizer)
0160   set_default_double_param("micromegas_enc", 670);
0161   set_default_double_param("micromegas_adc_threshold", 2680);
0162   set_default_double_param("micromegas_pedestal", 50000);
0163   set_default_double_param("micromegas_volts_per_charge", 20);
0164 }
0165
0166 //___________________________________________________________________________
0167 double PHG4MicromegasDigitizer::add_noise() const
0168 {
0169   return gsl_ran_gaussian(m_rng.get(), m_enc);
0170 }