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