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