ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHG4InttDigitizer.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHG4InttDigitizer.cc
1 // This is the new trackbase container version
2 
3 #include "PHG4InttDigitizer.h"
4 
5 #include "InttDeadMap.h"
6 
9 
10 // Move to new storage containers
11 #include <trackbase/TrkrDefs.h>
12 #include <trackbase/TrkrHit.h> // for TrkrHit
13 #include <trackbase/TrkrHitSet.h>
16 
17 #include <phparameter/PHParameterInterface.h> // for PHParameterInterface
18 
19 #include <intt/InttDefs.h>
20 
21 #include <fun4all/Fun4AllBase.h> // for Fun4AllBase::VERB...
23 #include <fun4all/SubsysReco.h> // for SubsysReco
24 
25 #include <phool/PHCompositeNode.h>
26 #include <phool/PHNode.h> // for PHNode
27 #include <phool/PHNodeIterator.h>
28 #include <phool/PHRandomSeed.h>
29 #include <phool/getClass.h>
30 #include <phool/phool.h> // for PHWHERE
31 
32 #include <TSystem.h>
33 
34 #include <gsl/gsl_randist.h>
35 #include <gsl/gsl_rng.h> // for gsl_rng_alloc
36 
37 #include <cassert>
38 #include <cfloat>
39 #include <cstdlib> // for exit
40 #include <iostream>
41 #include <memory> // for allocator_traits<...
42 #include <type_traits> // for __decay_and_strip...
43 
44 using namespace std;
45 
47  : SubsysReco(name)
48  , PHParameterInterface(name)
49  , detector("INTT")
50  , mNoiseMean(457.2)
51  , mNoiseSigma(166.6)
52  , mEnergyPerPair(3.62e-9) // GeV/e-h
53  , m_nCells(0)
54  , m_nDeadCells(0)
55 {
57  unsigned int seed = PHRandomSeed(); // fixed seed is handled in this funtcion
58  cout << Name() << " random seed: " << seed << endl;
59  RandomGenerator = gsl_rng_alloc(gsl_rng_mt19937);
60  gsl_rng_set(RandomGenerator, seed);
61 }
62 
64 {
65  gsl_rng_free(RandomGenerator);
66 }
67 
69 {
70  cout << "PHG4InttDigitizer::InitRun: detector = " << detector << endl;
71 
72  //-------------
73  // Add Hit Node
74  //-------------
75 
76  PHNodeIterator iter(topNode);
77 
78  // Looking for the DST node
79  PHCompositeNode *dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
80  if (!dstNode)
81  {
82  cout << PHWHERE << "DST Node missing, doing nothing." << endl;
84  }
85 
87 
88  // Create the run and par nodes
89  PHCompositeNode *runNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "RUN"));
90  PHCompositeNode *parNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "PAR"));
91 
92  string paramnodename = "G4CELLPARAM_" + detector;
93  string geonodename = "G4CELLGEO_" + detector;
94 
96  // save this to the run wise tree to store on DST
97  PHNodeIterator runIter(runNode);
98  PHCompositeNode *RunDetNode = dynamic_cast<PHCompositeNode *>(runIter.findFirst("PHCompositeNode", detector));
99  if (!RunDetNode)
100  {
101  RunDetNode = new PHCompositeNode(detector);
102  runNode->addNode(RunDetNode);
103  }
104  SaveToNodeTree(RunDetNode, paramnodename);
105  // save this to the parNode for use
106  PHNodeIterator parIter(parNode);
107  PHCompositeNode *ParDetNode = dynamic_cast<PHCompositeNode *>(parIter.findFirst("PHCompositeNode", detector));
108  if (!ParDetNode)
109  {
110  ParDetNode = new PHCompositeNode(detector);
111  parNode->addNode(ParDetNode);
112  }
113  PutOnParNode(ParDetNode, geonodename);
114  mNoiseMean = get_double_param("NoiseMean");
115  mNoiseSigma = get_double_param("NoiseSigma");
116  mEnergyPerPair = get_double_param("EnergyPerPair");
117 
118  //----------------
119  // Report Settings
120  //----------------
121 
122  if (Verbosity() > 0)
123  {
124  cout << "====================== PHG4InttDigitizer::InitRun() =====================" << endl;
125  for (std::map<int, unsigned int>::iterator iter1 = _max_adc.begin();
126  iter1 != _max_adc.end();
127  ++iter1)
128  {
129  cout << " Max ADC in Layer #" << iter1->first << " = " << iter1->second << endl;
130  }
131  for (std::map<int, float>::iterator iter2 = _energy_scale.begin();
132  iter2 != _energy_scale.end();
133  ++iter2)
134  {
135  cout << " Energy per ADC in Layer #" << iter2->first << " = " << 1.0e6 * iter2->second << " keV" << endl;
136  }
137  cout << "===========================================================================" << endl;
138  }
139 
141 }
142 
144 {
145  DigitizeLadderCells(topNode);
146 
148 }
149 
151 {
152  // FPHX 3-bit ADC, thresholds are set in "set_fphx_adc_scale".
153 
154  //PHG4CellContainer *cells = findNode::getClass<PHG4CellContainer>(topNode, "G4CELL_INTT");
155  PHG4CylinderGeomContainer *geom_container = findNode::getClass<PHG4CylinderGeomContainer>(topNode, "CYLINDERGEOM_INTT");
156 
157  //if (!geom_container || !cells) return;
158  if (!geom_container) return;
159 
160  PHG4CylinderGeomContainer::ConstRange layerrange = geom_container->get_begin_end();
161  for (PHG4CylinderGeomContainer::ConstIterator layeriter = layerrange.first;
162  layeriter != layerrange.second;
163  ++layeriter)
164  {
165  int layer = layeriter->second->get_layer();
166  if (_max_fphx_adc.find(layer) == _max_fphx_adc.end())
167  {
168  cout << "Error: _max_fphx_adc is not available." << endl;
169  gSystem->Exit(1);
170  }
171  float thickness = (layeriter->second)->get_thickness(); // cm
172  float mip_e = 0.003876 * thickness; // GeV
173  _energy_scale.insert(std::make_pair(layer, mip_e));
174  }
175 
176  return;
177 }
178 
180 {
181  //---------------------------
182  // Get common Nodes
183  //---------------------------
184  const InttDeadMap *deadmap = findNode::getClass<InttDeadMap>(topNode, "DEADMAP_INTT");
185  if (Verbosity() >= VERBOSITY_MORE)
186  {
187  if (deadmap)
188  {
189  cout << "PHG4InttDigitizer::DigitizeLadderCells - Use deadmap ";
190  deadmap->identify();
191  }
192  else
193  {
194  cout << "PHG4InttDigitizer::DigitizeLadderCells - Can not find deadmap, all channels enabled " << endl;
195  }
196  }
197 
198  // Get the TrkrHitSetContainer node
199  TrkrHitSetContainer *trkrhitsetcontainer = findNode::getClass<TrkrHitSetContainer>(topNode, "TRKR_HITSET");
200  if(!trkrhitsetcontainer)
201  {
202  cout << "Could not locate TRKR_HITSET node, quit! " << endl;
203  exit(1);
204  }
205 
206  // Get the TrkrHitTruthAssoc node
207  auto hittruthassoc = findNode::getClass<TrkrHitTruthAssoc>(topNode, "TRKR_HITTRUTHASSOC");
208 
209  //-------------
210  // Digitization
211  //-------------
212 
213  // We want all hitsets for the Intt
214  TrkrHitSetContainer::ConstRange hitset_range = trkrhitsetcontainer->getHitSets(TrkrDefs::TrkrId::inttId);
215  for (TrkrHitSetContainer::ConstIterator hitset_iter = hitset_range.first;
216  hitset_iter != hitset_range.second;
217  ++hitset_iter)
218  {
219  // we have an itrator to one TrkrHitSet for the intt from the trkrHitSetContainer
220  // get the hitset key so we can find the layer
221  TrkrDefs::hitsetkey hitsetkey = hitset_iter->first;
222  const int layer = TrkrDefs::getLayer(hitsetkey);
223  const int ladder_phi = InttDefs::getLadderPhiId(hitsetkey);
224  const int ladder_z = InttDefs::getLadderZId(hitsetkey);
225 
226  if(Verbosity() > 1)
227  cout << "PHG4InttDigitizer: found hitset with key: " << hitsetkey << " in layer " << layer << endl;
228 
229  // get all of the hits from this hitset
230  TrkrHitSet *hitset = hitset_iter->second;
231  TrkrHitSet::ConstRange hit_range = hitset->getHits();
232  std::set<TrkrDefs::hitkey> dead_hits; // hits on dead channel
233  for(TrkrHitSet::ConstIterator hit_iter = hit_range.first;
234  hit_iter != hit_range.second;
235  ++hit_iter)
236  {
237  ++m_nCells;
238 
239  TrkrHit *hit = hit_iter->second;
240  TrkrDefs::hitkey hitkey = hit_iter->first;
241  int strip_col = InttDefs::getCol(hitkey); // strip z index
242  int strip_row = InttDefs::getRow(hitkey); // strip phi index
243 
244  // Apply deadmap here if desired
245  if (deadmap)
246  {
247  if (deadmap->isDeadChannelIntt(
248  layer,
249  ladder_phi,
250  ladder_z,
251  strip_col,
252  strip_row
253  ))
254  {
255  ++m_nDeadCells;
256 
257  if (Verbosity() >= VERBOSITY_MORE)
258  {
259  cout << "PHG4InttDigitizer::DigitizeLadderCells - dead strip at layer " << layer << ": ";
260  hit->identify();
261  }
262 
263  dead_hits.insert(hit_iter->first); // store hitkey of dead channels to be remove later
264  continue;
265  }
266  } // if (deadmap)
267 
268  if (_energy_scale.count(layer) > 1)
269  assert(!"Error: _energy_scale has two or more keys.");
270 
271  const float mip_e = _energy_scale[layer];
272 
273  std::vector<std::pair<double, double> > vadcrange = _max_fphx_adc[layer];
274 
275  int adc = 0;
276  for (unsigned int irange = 0; irange < vadcrange.size(); ++irange)
277  if (hit->getEnergy() / TrkrDefs::InttEnergyScaleup >= vadcrange[irange].first * (double) mip_e && hit->getEnergy() / TrkrDefs::InttEnergyScaleup < vadcrange[irange].second * (double) mip_e)
278  adc = (unsigned short) irange;
279 
280  hit->setAdc(adc);
281 
282  if(Verbosity() > 2)
283  cout << "PHG4InttDigitizer: found hit with layer " << layer << " ladder_z " << ladder_z << " ladder_phi " << ladder_phi
284  << " strip_col " << strip_col << " strip_row " << strip_row << " adc " << hit->getAdc() << endl;
285 
286  } // end loop over hits in this hitset
287 
288  // remove hits on dead channel in TRKR_HITSET and TRKR_HITTRUTHASSOC
289  for(const auto& key:dead_hits) {
290  if(Verbosity() > 2)
291  cout<<" PHG4InttDigitizer: remove hit with key: " << key << endl;
292 
293  hitset->removeHit(key);
294 
295  if( hittruthassoc ) hittruthassoc->removeAssoc(hitsetkey, key);
296  }
297  } // end loop over hitsets
298 
299  return;
300 }
301 
304 {
305  if (Verbosity() >= VERBOSITY_SOME)
306  {
307  cout << "PHG4InttDigitizer::End - processed "
308  << m_nCells << " cell with "
309  << m_nDeadCells << " dead cells masked"
310  << " (" << 100. * m_nDeadCells / m_nCells << "%)" << endl;
311  }
312 
314 }
315 
317 {
318  set_default_double_param("NoiseMean", 457.2);
319  set_default_double_param("NoiseSigma", 166.6);
320  set_default_double_param("EnergyPerPair", 3.62e-9); // GeV/e-h
321  return;
322 }
323 
325 {
326 // float noise = gsl_ran_gaussian(RandomGenerator, mNoiseSigma) + mNoiseMean;
327 // noise = (noise < 0) ? 0 : noise;
328 
329  // Note the noise is bi-polar, i.e. can make ths signal fluctuate up and down.
330  // Much of the mNoiseSigma as extracted in https://github.com/sPHENIX-Collaboration/coresoftware/pull/580
331  // is statistical fluctuation from the limited calibration data. They does not directly apply here.
332  float noise = gsl_ran_gaussian(RandomGenerator, mNoiseMean);
333 
334  return noise;
335 }
336 
337 void PHG4InttDigitizer::set_adc_scale(const int &layer, const std::vector<double> &userrange)
338 {
339  if (userrange.size() != nadcbins)
340  {
341  cout << "Error: vector in set_fphx_adc_scale(vector) must have eight elements." << endl;
342  gSystem->Exit(1);
343  }
344  //sort(userrange.begin(), userrange.end()); // TODO, causes GLIBC error
345 
346  std::vector<std::pair<double, double> > vadcrange;
347  for (unsigned int irange = 0; irange < userrange.size(); ++irange)
348  {
349  if (irange == userrange.size() - 1)
350  {
351  vadcrange.push_back(std::make_pair(userrange[irange], FLT_MAX));
352  }
353  else
354  {
355  vadcrange.push_back(std::make_pair(userrange[irange], userrange[irange + 1]));
356  }
357  }
358  _max_fphx_adc.insert(std::make_pair(layer, vadcrange));
359 }