ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
RawTowerDigitizer.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file RawTowerDigitizer.cc
1 #include "RawTowerDigitizer.h"
2 
3 #include <calobase/RawTower.h> // for RawTower
4 #include <calobase/RawTowerContainer.h>
5 #include <calobase/RawTowerDeadMap.h>
6 #include <calobase/RawTowerDefs.h> // for keytype
7 #include <calobase/RawTowerGeom.h>
8 #include <calobase/RawTowerGeomContainer.h>
9 #include <calobase/RawTowerv1.h>
10 #include <calobase/RawTowerv2.h>
11 
12 #include <fun4all/Fun4AllBase.h> // for Fun4AllBase::VERBOSITY_MORE
14 #include <fun4all/SubsysReco.h> // for SubsysReco
15 
16 #include <phparameter/PHParameters.h>
17 
18 #include <phool/PHCompositeNode.h>
19 #include <phool/PHIODataNode.h>
20 #include <phool/PHNode.h> // for PHNode
21 #include <phool/PHNodeIterator.h>
22 #include <phool/PHObject.h> // for PHObject
23 #include <phool/PHRandomSeed.h>
24 #include <phool/getClass.h>
25 
26 #include <gsl/gsl_cdf.h>
27 #include <gsl/gsl_randist.h>
28 #include <gsl/gsl_rng.h>
29 
30 #include <cmath>
31 #include <cstdlib> // for exit
32 #include <exception> // for exception
33 #include <iostream>
34 #include <map>
35 #include <stdexcept>
36 #include <string>
37 #include <utility> // for pair
38 
39 using namespace std;
40 
42  : SubsysReco(name)
43  , m_DigiAlgorithm(kNo_digitization)
44  , m_SimTowers(nullptr)
45  , m_RawTowers(nullptr)
46  , m_RawTowerGeom(nullptr)
47  , m_DeadMap(nullptr)
48  , m_Detector("NONE")
49  , m_SimTowerNodePrefix("SIM")
50  , m_RawTowerNodePrefix("RAW")
51  , m_PhotonElecYieldVisibleGeV(NAN)
52  , m_PhotonElecADC(NAN)
53  , m_PedstalCentralADC(NAN)
54  , m_PedstalWidthADC(NAN)
55  , m_pedestalFile(false)
56  , m_ZeroSuppressionADC(-1000) //default to apply no zero suppression
57  , m_ZeroSuppressionFile(false)
58  , m_TowerType(-1)
59  , m_SiPMEffectivePixel(40000 * 4) // sPHENIX EMCal default, 4x Hamamatsu S12572-015P MPPC [sPHENIX TDR]
60  , _tower_params(name)
61 {
62  m_RandomGenerator = gsl_rng_alloc(gsl_rng_mt19937);
63  m_Seed = PHRandomSeed(); // fixed seed handled in PHRandomSeed()
64  // cout << Name() << " Random Seed: " << m_Seed << endl;
65  gsl_rng_set(m_RandomGenerator, m_Seed);
66 }
67 
69 {
70  gsl_rng_free(m_RandomGenerator);
71 }
72 
73 void RawTowerDigitizer::set_seed(const unsigned int iseed)
74 {
75  m_Seed = iseed;
76  gsl_rng_set(m_RandomGenerator, m_Seed);
77 }
78 
80 {
81  PHNodeIterator iter(topNode);
82 
83  // Looking for the DST node
84  PHCompositeNode *dstNode;
85  dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
86  if (!dstNode)
87  {
88  cout << Name() << "::" << m_Detector << "::" << __PRETTY_FUNCTION__
89  << "DST Node missing, doing nothing." << endl;
90  exit(1);
91  }
92 
93  try
94  {
95  CreateNodes(topNode);
96  }
97  catch (std::exception &e)
98  {
99  cout << e.what() << endl;
101  }
103 }
104 
106 {
107  if (Verbosity())
108  {
109  cout << Name() << "::" << m_Detector << "::" << __PRETTY_FUNCTION__
110  << "Process event entered. "
111  << "Digitalization method: ";
112 
114  {
115  cout << "directly pass the energy of sim tower to digitalized tower";
116  }
118  {
119  cout << "simple digitization with photon statistics, ADC conversion and pedstal";
120  }
121  cout << endl;
122  }
123  // loop over all possible towers, even empty ones. The digitization can add towers containing
124  // pedestals
126 
127  double deadChanEnergy = 0;
128 
129  for (RawTowerGeomContainer::ConstIterator it = all_towers.first;
130  it != all_towers.second; ++it)
131  {
132  const RawTowerDefs::keytype key = it->second->get_id();
134  const int eta = it->second->get_bineta();
135  const int phi = it->second->get_binphi();
136 
137  if (caloid == RawTowerDefs::LFHCAL)
138  {
139  const int l = it->second->get_binl();
140  if (m_ZeroSuppressionFile == true)
141  {
142  const string zsName = "ZS_ADC_eta" + to_string(eta) + "_phi" + to_string(phi) + "_l" + to_string(l);
145  }
146 
147  if (m_pedestalFile == true)
148  {
149  const string pedCentralName = "PedCentral_ADC_eta" + to_string(eta) + "_phi" + to_string(phi) + "_l" + to_string(l);
151  _tower_params.get_double_param(pedCentralName);
152  const string pedWidthName = "PedWidth_ADC_eta" + to_string(eta) + "_phi" + to_string(phi) + "_l" + to_string(l);
154  _tower_params.get_double_param(pedWidthName);
155  }
156  }
157  else
158  {
159  if (m_ZeroSuppressionFile == true)
160  {
161  const string zsName = "ZS_ADC_eta" + to_string(eta) + "_phi" + to_string(phi);
164  }
165 
166  if (m_pedestalFile == true)
167  {
168  const string pedCentralName = "PedCentral_ADC_eta" + to_string(eta) + "_phi" + to_string(phi);
170  _tower_params.get_double_param(pedCentralName);
171  const string pedWidthName = "PedWidth_ADC_eta" + to_string(eta) + "_phi" + to_string(phi);
173  _tower_params.get_double_param(pedWidthName);
174  }
175  }
176 
177 
178 
179  if (m_TowerType >= 0)
180  {
181  // Skip towers that don't match the type we are supposed to digitize
182  if (m_TowerType != it->second->get_tower_type())
183  {
184  continue;
185  }
186  }
187 
188  RawTower *sim_tower = m_SimTowers->getTower(key);
189  if (m_DeadMap)
190  {
191  if (m_DeadMap->isDeadTower(key))
192  {
193  if (sim_tower) deadChanEnergy += sim_tower->get_energy();
194 
195  sim_tower = nullptr;
196 
197  if (Verbosity() >= VERBOSITY_MORE)
198  {
199  cout << Name() << "::" << m_Detector << "::" << __PRETTY_FUNCTION__
200  << " apply dead tower " << key << endl;
201  }
202  }
203  }
204 
205  RawTower *digi_tower = nullptr;
206 
208  {
209  // for no digitization just copy existing towers
210  if (sim_tower)
211  {
212  digi_tower = new RawTowerv2(*sim_tower);
213  }
214  }
216  {
217  // for photon digitization towers can be created if sim_tower is null pointer
218  digi_tower = simple_photon_digitization(sim_tower);
219  }
221  {
222  // for photon digitization towers can be created if sim_tower is null pointer
223  digi_tower = sipm_photon_digitization(sim_tower);
224  }
225  else
226  {
227  cout << Name() << "::" << m_Detector << "::" << __PRETTY_FUNCTION__
228  << " invalid digitization algorithm #" << m_DigiAlgorithm
229  << endl;
230 
232  }
233 
234  if (digi_tower)
235  {
236  m_RawTowers->AddTower(key, digi_tower);
237 
238  if (Verbosity() >= VERBOSITY_MORE)
239  {
240  cout << Name() << "::" << m_Detector << "::" << __PRETTY_FUNCTION__
241  << " output tower:"
242  << endl;
243  digi_tower->identify();
244  }
245  }
246  }
247 
248  if (Verbosity())
249  {
250  cout << Name() << "::" << m_Detector << "::" << __PRETTY_FUNCTION__
251  << "input sum energy = " << m_SimTowers->getTotalEdep() << " GeV"
252  << ", dead channel masked energy = " << deadChanEnergy << " GeV"
253  << ", output sum digitalized value = " << m_RawTowers->getTotalEdep() << " ADC"
254  << endl;
255  }
257 }
258 
259 RawTower *
261 {
262  RawTower *digi_tower = nullptr;
263 
264  double energy = 0;
265  if (sim_tower)
266  {
267  energy = sim_tower->get_energy();
268  }
269  const double photon_count_mean = energy * m_PhotonElecYieldVisibleGeV;
270  const int photon_count = gsl_ran_poisson(m_RandomGenerator, photon_count_mean);
271  const int signal_ADC = floor(photon_count / m_PhotonElecADC);
272 
273  const double pedstal = m_PedstalCentralADC + ((m_PedstalWidthADC > 0) ? gsl_ran_gaussian(m_RandomGenerator, m_PedstalWidthADC) : 0);
274  const int sum_ADC = signal_ADC + (int) pedstal;
275 
276  if (sum_ADC > m_ZeroSuppressionADC)
277  {
278  // create new digitalizaed tower
279  if (sim_tower)
280  {
281  digi_tower = new RawTowerv2(*sim_tower);
282  }
283  else
284  {
285  digi_tower = new RawTowerv2();
286  }
287  digi_tower->set_energy((double) sum_ADC);
288  }
289 
290  if (Verbosity() >= 2)
291  {
292  cout << Name() << "::" << m_Detector << "::" << __PRETTY_FUNCTION__
293  << endl;
294 
295  cout << "input: ";
296  if (sim_tower)
297  {
298  sim_tower->identify();
299  }
300  else
301  {
302  cout << "None" << endl;
303  }
304  cout << "output based on "
305  << "sum_ADC = " << sum_ADC << ", zero_sup = "
306  << m_ZeroSuppressionADC << " : ";
307  if (digi_tower)
308  {
309  digi_tower->identify();
310  }
311  else
312  {
313  cout << "None" << endl;
314  }
315  }
316 
317  return digi_tower;
318 }
319 
320 RawTower *
322 {
323  RawTower *digi_tower = nullptr;
324 
325  double energy = 0;
326  if (sim_tower)
327  {
328  energy = sim_tower->get_energy();
329  }
330 
331  int signal_ADC = 0;
332 
333  if (energy > 0)
334  {
335  const double photon_count_mean = energy * m_PhotonElecYieldVisibleGeV;
336  const double poission_param_per_pixel = photon_count_mean / m_SiPMEffectivePixel;
337  const double prob_activated_per_pixel = gsl_cdf_poisson_Q(0, poission_param_per_pixel);
338  const double active_pixel = gsl_ran_binomial(m_RandomGenerator, prob_activated_per_pixel, m_SiPMEffectivePixel);
339  signal_ADC = floor(active_pixel / m_PhotonElecADC);
340  }
341 
342  const double pedstal = m_PedstalCentralADC + ((m_PedstalWidthADC > 0) ? gsl_ran_gaussian(m_RandomGenerator, m_PedstalWidthADC) : 0);
343  const int sum_ADC = signal_ADC + (int) pedstal;
344 
345  if (sum_ADC > m_ZeroSuppressionADC)
346  {
347  // create new digitalizaed tower
348  if (sim_tower)
349  {
350  digi_tower = new RawTowerv2(*sim_tower);
351  }
352  else
353  {
354  digi_tower = new RawTowerv2();
355  }
356  digi_tower->set_energy((double) sum_ADC);
357  }
358 
359  if (Verbosity() >= 2)
360  {
361  cout << Name() << "::" << m_Detector << "::" << __PRETTY_FUNCTION__
362  << endl;
363 
364  cout << "input: ";
365  if (sim_tower)
366  {
367  sim_tower->identify();
368  }
369  else
370  {
371  cout << "None" << endl;
372  }
373  cout << "output based on "
374  << "sum_ADC = " << sum_ADC << ", zero_sup = "
375  << m_ZeroSuppressionADC << " : ";
376  if (digi_tower)
377  {
378  digi_tower->identify();
379  }
380  else
381  {
382  cout << "None" << endl;
383  }
384  }
385 
386  return digi_tower;
387 }
388 
390 {
391  PHNodeIterator iter(topNode);
392  PHCompositeNode *runNode = static_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "RUN"));
393  if (!runNode)
394  {
395  cout << Name() << "::" << m_Detector << "::" << __PRETTY_FUNCTION__
396  << "Run Node missing, doing nothing." << endl;
397  throw std::runtime_error(
398  "Failed to find Run node in RawTowerDigitizer::CreateNodes");
399  }
400 
401  string TowerGeomNodeName = "TOWERGEOM_" + m_Detector;
402  m_RawTowerGeom = findNode::getClass<RawTowerGeomContainer>(topNode, TowerGeomNodeName);
403  if (!m_RawTowerGeom)
404  {
405  cout << Name() << "::" << m_Detector << "::" << __PRETTY_FUNCTION__
406  << " " << TowerGeomNodeName << " Node missing, doing bail out!"
407  << endl;
408  throw std::runtime_error("Failed to find " + TowerGeomNodeName + " node in RawTowerDigitizer::CreateNodes");
409  }
410 
411  if (Verbosity() >= 1)
412  {
414  }
415 
416  const string deadMapName = "DEADMAP_" + m_Detector;
417  m_DeadMap = findNode::getClass<RawTowerDeadMap>(topNode, deadMapName);
418  if (m_DeadMap)
419  {
420  cout << Name() << "::" << m_Detector << "::" << __PRETTY_FUNCTION__
421  << " use dead map: ";
422  m_DeadMap->identify();
423  }
424 
425  PHCompositeNode *dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
426  if (!dstNode)
427  {
428  cout << Name() << "::" << m_Detector << "::" << __PRETTY_FUNCTION__
429  << "DST Node missing, doing nothing." << endl;
430  throw std::runtime_error("Failed to find DST node in RawTowerDigitizer::CreateNodes");
431  }
432 
433  string SimTowerNodeName = "TOWER_" + m_SimTowerNodePrefix + "_" + m_Detector;
434  m_SimTowers = findNode::getClass<RawTowerContainer>(dstNode, SimTowerNodeName);
435  if (!m_SimTowers)
436  {
437  cout << Name() << "::" << m_Detector << "::" << __PRETTY_FUNCTION__
438  << " " << SimTowerNodeName << " Node missing, doing bail out!"
439  << endl;
440  throw std::runtime_error("Failed to find " + SimTowerNodeName + " node in RawTowerDigitizer::CreateNodes");
441  }
442 
443  // Create the tower nodes on the tree
444  PHNodeIterator dstiter(dstNode);
445  PHCompositeNode *DetNode = dynamic_cast<PHCompositeNode *>(dstiter.findFirst("PHCompositeNode", m_Detector));
446  if (!DetNode)
447  {
448  DetNode = new PHCompositeNode(m_Detector);
449  dstNode->addNode(DetNode);
450  }
451 
452  // Be careful as a previous digitizer may have been registered for this detector
453  string RawTowerNodeName = "TOWER_" + m_RawTowerNodePrefix + "_" + m_Detector;
454  m_RawTowers = findNode::getClass<RawTowerContainer>(DetNode, RawTowerNodeName);
455  if (!m_RawTowers)
456  {
459  RawTowerNodeName, "PHObject");
460  DetNode->addNode(towerNode);
461  }
462 
463  return;
464 }