Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Prototype3DSTReader.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file Prototype3DSTReader.cc
1 // $Id: Prototype3DSTReader.cc,v 1.11 2015/01/06 02:52:07 jinhuang Exp $
2 
11 #include "Prototype3DSTReader.h"
12 
13 #include <calobase/RawTower.h> // for RawTower
14 #include <calobase/RawTowerContainer.h>
15 
16 #include <pdbcalbase/PdbParameterMap.h>
17 
18 #include <phparameter/PHParameters.h>
19 
20 #include <fun4all/PHTFileServer.h>
21 #include <fun4all/SubsysReco.h> // for SubsysReco
22 
23 #include <phool/getClass.h>
24 
25 #include <TClass.h> // for TClass
26 #include <TClonesArray.h> // for TClonesArray
27 #include <TObject.h> // for TObject
28 #include <TString.h> // for Form
29 #include <TTree.h>
30 
31 #include <cassert>
32 #include <iostream> // for operator<<, basic_ostream, endl
33 #include <map>
34 #include <utility> // for pair
35 
36 class PHCompositeNode;
37 
38 using namespace std;
39 
41  : SubsysReco("Prototype3DSTReader")
42  , nblocks(0)
43  , _event(0)
44  , //
45  _out_file_name(filename)
46  , /*_file(nullptr), */ _T(nullptr)
47  , //
48  _tower_zero_sup(-10000000)
49 {
50 }
51 
53 {
54  cout << "Prototype3DSTReader::destructor - Clean ups" << endl;
55 
56  if (_T)
57  {
58  _T->ResetBranchAddresses();
59  }
60 
61  _records.clear();
62 }
63 
65 {
66  const static int arr_size = 100;
67 
68  if (_tower_postfix.size())
69  {
70  cout << "Prototype3DSTReader::Init - zero suppression for calorimeter "
71  "towers = "
72  << _tower_zero_sup << " GeV" << endl;
73  }
74  for (vector<string>::const_iterator it = _runinfo_list.begin();
75  it != _runinfo_list.end(); ++it)
76  {
77  const string &nodenam = *it;
78 
79  record rec;
80  rec._cnt = 0;
81  rec._name = nodenam;
82  rec._arr = nullptr;
83  rec._arr_ptr = nullptr;
84  rec._dvalue = 0;
86 
87  _records.push_back(rec);
88 
89  nblocks++;
90  }
91  for (vector<string>::const_iterator it = _eventinfo_list.begin();
92  it != _eventinfo_list.end(); ++it)
93  {
94  const string &nodenam = *it;
95 
96  record rec;
97  rec._cnt = 0;
98  rec._name = nodenam;
99  rec._arr = nullptr;
100  rec._arr_ptr = nullptr;
101  rec._dvalue = 0;
103 
104  _records.push_back(rec);
105 
106  nblocks++;
107  }
108 
109  for (vector<string>::const_iterator it = _tower_postfix.begin();
110  it != _tower_postfix.end(); ++it)
111  {
112  const char *class_name = RawTower_type::Class()->GetName();
113 
114  const string &nodenam = *it;
115 
116  string hname = Form("TOWER_%s", nodenam.c_str());
117  // _node_name.push_back(hname);
118  cout << "Prototype3DSTReader::Init - saving raw tower info from node: "
119  << hname << " - " << class_name << endl;
120 
121  record rec;
122  rec._cnt = 0;
123  rec._name = hname;
124  rec._arr = make_shared<TClonesArray>(class_name, arr_size);
125  rec._arr_ptr = rec._arr.get();
126  rec._dvalue = 0;
127  rec._type = record::typ_tower;
128 
129  _records.push_back(rec);
130 
131  nblocks++;
132  }
133 
134  for (vector<string>::const_iterator it = _towertemp_postfix.begin();
135  it != _towertemp_postfix.end(); ++it)
136  {
137  const string &nodenam = *it;
138  string hname = Form("TOWER_TEMPERATURE_%s", nodenam.c_str());
139 
140  cout << "Prototype3DSTReader::Init - saving average tower temperature info "
141  "from node: "
142  << hname << endl;
143 
144  record rec;
145  rec._cnt = 0;
146  rec._name = hname;
147  rec._arr = nullptr;
148  rec._arr_ptr = nullptr;
149  rec._dvalue = 0;
151 
152  _records.push_back(rec);
153 
154  nblocks++;
155  }
156  cout << "Prototype3DSTReader::Init - requested " << nblocks << " nodes"
157  << endl;
158 
159  build_tree();
160 
161  return 0;
162 }
163 
165 {
166  cout << "Prototype3DSTReader::build_tree - output to " << _out_file_name
167  << endl;
168 
169  static const int BUFFER_SIZE = 32000;
170 
171  // open TFile
172  PHTFileServer::get().open(_out_file_name, "RECREATE");
173 
174  _T = new TTree("T", "Prototype3DSTReader");
175 
176  nblocks = 0;
177  for (records_t::iterator it = _records.begin(); it != _records.end(); ++it)
178  {
179  record &rec = *it;
180 
181  cout << "Prototype3DSTReader::build_tree - Add " << rec._name << endl;
182 
183  if (rec._type == record::typ_runinfo)
184  {
185  const string name_cnt = rec._name;
186  const string name_cnt_desc = name_cnt + "/D";
187  _T->Branch(name_cnt.c_str(), &(rec._dvalue), name_cnt_desc.c_str(),
188  BUFFER_SIZE);
189  }
190  if (rec._type == record::typ_eventinfo)
191  {
192  const string name_cnt = rec._name;
193  const string name_cnt_desc = name_cnt + "/D";
194  _T->Branch(name_cnt.c_str(), &(rec._dvalue), name_cnt_desc.c_str(),
195  BUFFER_SIZE);
196  }
197  else if (rec._type == record::typ_tower)
198  {
199  const string name_cnt = "n_" + rec._name;
200  const string name_cnt_desc = name_cnt + "/I";
201  _T->Branch(name_cnt.c_str(), &(rec._cnt), name_cnt_desc.c_str(),
202  BUFFER_SIZE);
203  _T->Branch(rec._name.c_str(), &(rec._arr_ptr), BUFFER_SIZE, 99);
204  }
205  else if (rec._type == record::typ_towertemp)
206  {
207  const string name_cnt = rec._name + "_AVG";
208  const string name_cnt_desc = name_cnt + "/D";
209  _T->Branch(name_cnt.c_str(), &(rec._dvalue), name_cnt_desc.c_str(),
210  BUFFER_SIZE);
211  }
212 
213  nblocks++;
214  }
215 
216  cout << "Prototype3DSTReader::build_tree - added " << nblocks << " nodes"
217  << endl;
218 
219  _T->SetAutoSave(16000);
220 }
221 
223 {
224  // const double significand = _event / TMath::Power(10, (int)
225  // (log10(_event)));
226  //
227  // if (fmod(significand, 1.0) == 0 && significand <= 10)
228  // cout << "Prototype3DSTReader::process_event - " << _event << endl;
229  _event++;
230 
231  for (records_t::iterator it = _records.begin(); it != _records.end(); ++it)
232  {
233  record &rec = *it;
234 
235  rec._cnt = 0;
236 
237  if (rec._type == record::typ_hit)
238  {
239  assert(0);
240  } // if (rec._type == record::typ_hit)
241  else if (rec._type == record::typ_tower)
242  {
243  assert(rec._arr.get() == rec._arr_ptr);
244  assert(rec._arr.get());
245  rec._arr->Clear();
246 
247  if (Verbosity() >= 2)
248  cout << "Prototype3DSTReader::process_event - processing tower "
249  << rec._name << endl;
250 
251  RawTowerContainer *hits =
252  findNode::getClass<RawTowerContainer>(topNode, rec._name);
253  if (!hits)
254  {
255  if (_event < 2)
256  cout << "Prototype3DSTReader::process_event - Error - can not find "
257  "node "
258  << rec._name << endl;
259  }
260  else
261  {
262  RawTowerContainer::ConstRange hit_range = hits->getTowers();
263 
264  if (Verbosity() >= 2)
265  cout << "Prototype3DSTReader::process_event - processing "
266  << rec._name << " and received " << hits->size() << " tower hits"
267  << endl;
268 
269  for (RawTowerContainer::ConstIterator hit_iter = hit_range.first;
270  hit_iter != hit_range.second; hit_iter++)
271  {
272  RawTower *hit_raw = hit_iter->second;
273 
274  RawTower_type *hit = dynamic_cast<RawTower_type *>(hit_raw);
275  // RawTower * hit = hit_iter->second;
276 
277  assert(hit);
278 
279  if (hit->get_energy() < _tower_zero_sup)
280  {
281  if (Verbosity() >= 2)
282  cout << "Prototype3DSTReader::process_event - suppress low "
283  "energy tower hit "
284  << rec._name
285  << " @ ("
286  // << hit->get_thetaMin()
287  // << ", " << hit->get_phiMin()
288  << "), Energy = " << hit->get_energy() << endl;
289 
290  continue;
291  }
292 
293  new ((*(rec._arr.get()))[rec._cnt]) RawTower_type();
294 
295  if (Verbosity() >= 2)
296  cout << "Prototype3DSTReader::process_event - processing Tower hit "
297  << rec._name
298  << " @ ("
299  // << hit->get_thetaMin() << ", "
300  // << hit->get_phiMin()
301  << "), Energy = " << hit->get_energy() << " - "
302  << rec._arr.get()->At(rec._cnt)->ClassName() << endl;
303 
304  // rec._arr->Print();
305 
306  RawTower_type *new_hit =
307  dynamic_cast<RawTower_type *>(rec._arr.get()->At(rec._cnt));
308  assert(new_hit);
309 
310  *new_hit = (*hit);
311 
312  rec._cnt++;
313  }
314  } // if (!hits)
315  } // if (rec._type == record::typ_hit)
316  else if (rec._type == record::typ_towertemp)
317  {
318  if (Verbosity() >= 2)
319  cout << "Prototype3DSTReader::process_event - processing tower "
320  "temperature "
321  << rec._name << endl;
322 
323  RawTowerContainer *hits =
324  findNode::getClass<RawTowerContainer>(topNode, rec._name);
325  if (!hits)
326  {
327  if (_event < 2)
328  cout << "Prototype3DSTReader::process_event - Error - can not find "
329  "node "
330  << rec._name << endl;
331  }
332  else
333  {
334  RawTowerContainer::ConstRange hit_range = hits->getTowers();
335 
336  if (Verbosity() >= 2)
337  cout << "Prototype3DSTReader::process_event - processing "
338  << rec._name << " and received " << hits->size() << " tower hits"
339  << endl;
340 
341  rec._cnt = 0;
342 
343  for (RawTowerContainer::ConstIterator hit_iter = hit_range.first;
344  hit_iter != hit_range.second; hit_iter++)
345  {
346  RawTower *hit_raw = hit_iter->second;
347 
348  RawTowerT_type *hit = dynamic_cast<RawTowerT_type *>(hit_raw);
349  // RawTower * hit = hit_iter->second;
350 
351  assert(hit);
352 
353  rec._dvalue +=
355 
356  ++rec._cnt;
357  }
358 
359  rec._dvalue /= rec._cnt;
360  } // if (!hits)
361  } // if (rec._type == record::typ_hit)
362  else if (rec._type == record::typ_runinfo)
363  {
364  PdbParameterMap *info =
365  findNode::getClass<PdbParameterMap>(topNode, "RUN_INFO");
366 
367  assert(info);
368 
369  PHParameters run_info_copy("RunInfo");
370  run_info_copy.FillFrom(info);
371 
372  rec._dvalue = run_info_copy.get_double_param(rec._name);
373 
374  } //
375  else if (rec._type == record::typ_eventinfo)
376  {
377  PdbParameterMap *info =
378  findNode::getClass<PdbParameterMap>(topNode, "EVENT_INFO");
379 
380  assert(info);
381 
382  PHParameters event_info_copy("EVENT_INFO");
383  event_info_copy.FillFrom(info);
384 
385  rec._dvalue = event_info_copy.get_double_param(rec._name);
386 
387  } // if (rec._type == record::typ_hit)
388  else if (rec._type == record::typ_part)
389  {
390  assert(0);
391  } // if (rec._type == record::typ_part)
392  else if (rec._type == record::typ_vertex)
393  {
394  assert(0);
395  } // if (_load_all_particle)
396 
397  } // for (records_t::iterator it = _records.begin(); it != _records.end();
398  // ++it)
399 
400  if (_T)
401  _T->Fill();
402 
403  return 0;
404 } // for (records_t::iterator it = _records.begin(); it != _records.end();
405 // ++it)
406 
408 {
409  cout << "Prototype3DSTReader::End - Clean ups" << endl;
410 
411  if (_T)
412  {
414  _T->Write();
415  _T->ResetBranchAddresses();
416  }
417 
418  _records.clear();
419 
420  return 0;
421 }