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