ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHG4MicromegasDigitizer.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHG4MicromegasDigitizer.cc
1 
7 
8 // Move to new storage containers
9 #include <trackbase/TrkrDefs.h>
10 #include <trackbase/TrkrHit.h>
11 #include <trackbase/TrkrHitSet.h>
14 
15 #include <phparameter/PHParameterInterface.h> // for PHParameterInterface
16 
18 #include <fun4all/SubsysReco.h>
19 
20 #include <phool/PHRandomSeed.h>
21 #include <phool/getClass.h>
22 
23 #include <gsl/gsl_randist.h>
24 #include <gsl/gsl_rng.h> // for gsl_rng_alloc, gsl_rng...
25 
26 #include <cassert>
27 #include <iostream> // for operator<<, basic_ostream
28 #include <set>
29 #include <utility> // for pair
30 
31 namespace
32 {
33 
34  // local version of std::clamp, which is only available for c++17
35  template<class T>
36  constexpr const T& clamp( const T& v, const T& lo, const T& hi )
37  { return (v < lo) ? lo : (hi < v) ? hi : v; }
38 
39 }
40 //____________________________________________________________________________
42  : SubsysReco(name)
43  , PHParameterInterface(name)
44 {
45  // initialize rng
46  const unsigned int seed = PHRandomSeed();
47  m_rng.reset( gsl_rng_alloc(gsl_rng_mt19937) );
48  gsl_rng_set( m_rng.get(), seed );
49 
51 }
52 
53 //____________________________________________________________________________
55 {
56 
58 
59  // load parameters
60  m_adc_threshold = get_double_param( "micromegas_adc_threshold" );
61  m_enc = get_double_param( "micromegas_enc" );
62  m_pedestal = get_double_param( "micromegas_pedestal" );
63  m_volts_per_charge = get_double_param( "micromegas_volts_per_charge" );
64 
65  // printout
66  std::cout
67  << "PHG4MicromegasDigitizer::InitRun\n"
68  << " m_adc_threshold: " << m_adc_threshold << " electrons\n"
69  << " m_enc: " << m_enc << " electrons\n"
70  << " m_pedestal: " << m_pedestal << " electrons\n"
71  << " m_volts_per_charge: " << m_volts_per_charge << " mV/fC\n"
72  << std::endl;
73 
74  // threshold is effectively applied on top of pedestal
76 
77  /*
78  * Factor that convertes charge in a voltage in each z bin
79  * the scale up factor of 2.4 is meant to account for shaping time (80ns)
80  * it only applies to the signal
81  * see: simulations/g4simulations/g4tpc/PHG4TpcDigitizer::InitRun
82  */
85 
87 }
88 
89 //____________________________________________________________________________
91 {
92  // Get Nodes
93  // Get the TrkrHitSetContainer node
94  auto trkrhitsetcontainer = findNode::getClass<TrkrHitSetContainer>(topNode, "TRKR_HITSET");
95  assert( trkrhitsetcontainer );
96 
97  // Get the TrkrHitTruthAssoc node
98  auto hittruthassoc = findNode::getClass<TrkrHitTruthAssoc>(topNode, "TRKR_HITTRUTHASSOC");
99 
100  // get all micromegas hitsets
101  const auto hitset_range = trkrhitsetcontainer->getHitSets(TrkrDefs::TrkrId::micromegasId);
102  for( auto hitset_it = hitset_range.first; hitset_it != hitset_range.second; ++hitset_it )
103  {
104 
105  // get key
106  const TrkrDefs::hitsetkey hitsetkey = hitset_it->first;
107 
108  // get all of the hits from this hitset
109  TrkrHitSet* hitset = hitset_it->second;
110  TrkrHitSet::ConstRange hit_range = hitset->getHits();
111 
112  // keep track of hits to be removed
113  std::set<TrkrDefs::hitkey> removed_keys;
114 
115  // loop over hits
116  for( auto hit_it = hit_range.first; hit_it != hit_range.second; ++hit_it )
117  {
118  // store key and hit
119  const TrkrDefs::hitkey& key = hit_it->first;
120  TrkrHit *hit = hit_it->second;
121 
122  // get energy (electrons)
123  const double signal = hit->getEnergy();
124  const double noise = add_noise();
125 
126  // convert to mV
127  const double voltage = (m_pedestal + noise)*m_volt_per_electron_noise + signal*m_volt_per_electron_signal;
128 
129  // compare to threshold
131  {
132  // keep hit, update adc
133  hit->setAdc( clamp<uint>( voltage*m_adc_per_volt, 0, 1023 ) );
134  } else {
135  // mark hit as removable
136  removed_keys.insert( key );
137  }
138 
139  }
140 
141  // remove hits
142  for( const auto& key:removed_keys )
143  {
144  hitset->removeHit(key);
145  if( hittruthassoc ) hittruthassoc->removeAssoc(hitsetkey, key);
146  }
147 
148  }
150 }
151 
152 //___________________________________________________________________________
154 {
155  // all values taken from TPC sampa chips (simulations/g4simulations/g4tpc/PHG4TpcDigitizer)
156  set_default_double_param( "micromegas_enc", 670 );
157  set_default_double_param( "micromegas_adc_threshold", 2680 );
158  set_default_double_param( "micromegas_pedestal", 50000 );
159  set_default_double_param( "micromegas_volts_per_charge", 20 );
160 }
161 
162 //___________________________________________________________________________
164 { return gsl_ran_gaussian( m_rng.get(), m_enc); }