EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHG4DSTReader.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHG4DSTReader.cc
1 // $Id: PHG4DSTReader.cc,v 1.11 2015/01/06 02:52:07 jinhuang Exp $
2 
11 #include "PHG4DSTReader.h"
12 
13 #include <g4main/PHG4HitContainer.h>
14 #include <g4main/PHG4TruthInfoContainer.h>
15 
16 #include <calobase/RawTower.h>
17 #include <calobase/RawTowerContainer.h>
18 #include <calobase/RawTowerv1.h>
19 
20 #include <g4main/PHG4Hit.h>
21 #include <g4main/PHG4HitEval.h>
22 #include <g4main/PHG4Particle.h>
23 #include <g4main/PHG4Particlev2.h>
24 #include <g4main/PHG4VtxPoint.h>
25 #include <g4main/PHG4VtxPointv1.h>
26 
27 #include <g4jets/Jet.h>
28 #include <g4jets/JetMap.h>
29 #include <g4jets/Jetv1.h>
30 
31 #include <fun4all/Fun4AllReturnCodes.h>
32 #include <fun4all/PHTFileServer.h>
33 #include <fun4all/SubsysReco.h>
34 
35 #include <phool/getClass.h>
36 
37 #include <TClass.h>
38 #include <TClonesArray.h>
39 #include <TObject.h>
40 #include <TString.h>
41 #include <TTree.h>
42 
43 #include <cassert>
44 #include <iostream>
45 #include <set>
46 #include <sstream>
47 #include <utility>
48 
49 using namespace std;
50 
51 typedef PHG4Particlev2 part_type;
52 typedef PHG4HitEval hit_type;
53 typedef PHG4VtxPointv1 vertex_type;
54 typedef RawTowerv1 RawTower_type;
55 typedef Jetv1 PHPyJet_type;
56 
57 PHG4DSTReader::PHG4DSTReader(const string &filename)
58  : SubsysReco("PHG4DSTReader")
59  , nblocks(0)
60  , _event(0)
61  , //
62  _out_file_name(filename)
63  , /*_file(nullptr), */ _T(nullptr)
64  , //
65  _save_particle(true)
66  , _load_all_particle(false)
67  , _load_active_particle(
68  true)
69  , _save_vertex(true)
70  , _tower_zero_sup(.0)
71 {
72  // TODO Auto-generated constructor stub
73 }
74 
75 PHG4DSTReader::~PHG4DSTReader()
76 {
77  if (Verbosity() > 1)
78  {
79  cout << "PHG4DSTReader::destructor - Clean ups" << endl;
80  }
81  if (_T)
82  {
83  _T->ResetBranchAddresses();
84  }
85 
86  _records.clear();
87 }
88 
89 int PHG4DSTReader::Init(PHCompositeNode *)
90 {
91  const int arr_size = 100;
92 
93  // BOOST_FOREACH(string nodenam, _node_postfix)
94  for (vector<string>::const_iterator it = _node_postfix.begin();
95  it != _node_postfix.end(); ++it)
96  {
97  const char *class_name = hit_type::Class()->GetName();
98 
99  const string &nodenam = *it;
100 
101  string hname = Form("G4HIT_%s", nodenam.c_str());
102  // _node_name.push_back(hname);
103  if (Verbosity() > 0)
104  {
105  cout << "PHG4DSTReader::Init - saving hits from node: " << hname << " - "
106  << class_name << endl;
107  }
108 
109  record rec;
110  rec._cnt = 0;
111  rec._name = hname;
112  rec._arr = boost::make_shared<TClonesArray>(class_name, arr_size);
113  rec._arr_ptr = rec._arr.get();
114  rec._type = record::typ_hit;
115 
116  _records.push_back(rec);
117 
118  nblocks++;
119  }
120 
121  if (_tower_postfix.size() && Verbosity() > 0)
122  {
123  cout << "PHG4DSTReader::Init - zero suppression for calorimeter towers = "
124  << _tower_zero_sup << " GeV" << endl;
125  }
126  for (vector<string>::const_iterator it = _tower_postfix.begin();
127  it != _tower_postfix.end(); ++it)
128  {
129  const char *class_name = RawTower_type::Class()->GetName();
130 
131  const string &nodenam = *it;
132 
133  string hname = Form("TOWER_%s", nodenam.c_str());
134  // _node_name.push_back(hname);
135  if (Verbosity() > 0)
136  {
137  cout << "PHG4DSTReader::Init - saving raw tower info from node: " << hname
138  << " - " << class_name << endl;
139  }
140  record rec;
141  rec._cnt = 0;
142  rec._name = hname;
143  rec._arr = boost::make_shared<TClonesArray>(class_name, arr_size);
144  rec._arr_ptr = rec._arr.get();
145  rec._type = record::typ_tower;
146 
147  _records.push_back(rec);
148 
149  nblocks++;
150  }
151 
152  for (vector<string>::const_iterator it = _jet_postfix.begin();
153  it != _jet_postfix.end(); ++it)
154  {
155  const char *class_name = PHPyJet_type::Class()->GetName();
156 
157  const string &hname = *it;
158  if (Verbosity() > 0)
159  {
160  cout << "PHG4DSTReader::Init - saving jets from node: " << hname << " - "
161  << class_name << endl;
162  }
163  record rec;
164  rec._cnt = 0;
165  rec._name = hname;
166  rec._arr = boost::make_shared<TClonesArray>(class_name, arr_size);
167  rec._arr_ptr = rec._arr.get();
168  rec._type = record::typ_jets;
169 
170  _records.push_back(rec);
171 
172  nblocks++;
173  }
174 
175  if (_save_particle)
176  {
177  //save particles
178 
179  const char *class_name = part_type::Class()->GetName();
180 
181  if (Verbosity() > 0)
182  {
183  cout << "PHG4DSTReader::Init - saving Particles node:" << class_name
184  << endl;
185  }
186  record rec;
187  rec._cnt = 0;
188  rec._name = "PHG4Particle";
189  rec._arr = boost::make_shared<TClonesArray>(class_name, arr_size);
190  rec._arr_ptr = rec._arr.get();
191  rec._type = record::typ_part;
192 
193  _records.push_back(rec);
194 
195  nblocks++;
196  }
197 
198  if (_save_vertex)
199  {
200  //save particles
201 
202  const char *class_name = vertex_type::Class()->GetName();
203 
204  if (Verbosity() > 0)
205  {
206  cout << "PHG4DSTReader::Init - saving vertex for particles under node:"
207  << class_name << endl;
208  }
209  record rec;
210  rec._cnt = 0;
211  rec._name = "PHG4VtxPoint";
212  rec._arr = boost::make_shared<TClonesArray>(class_name, arr_size);
213  rec._arr_ptr = rec._arr.get();
214  rec._type = record::typ_vertex;
215 
216  _records.push_back(rec);
217 
218  nblocks++;
219  }
220  if (Verbosity() > 0)
221  {
222  cout << "PHG4DSTReader::Init - requested " << nblocks << " nodes" << endl;
223  }
224  build_tree();
225 
226  return 0;
227 }
228 
229 void PHG4DSTReader::build_tree()
230 {
231  if (Verbosity() > 0)
232  {
233  cout << "PHG4DSTReader::build_tree - output to " << _out_file_name << endl;
234  }
235  static const int BUFFER_SIZE = 32000;
236 
237  // open TFile
238  PHTFileServer::get().open(_out_file_name, "RECREATE");
239 
240  _T = new TTree("T", "PHG4DSTReader");
241 
242  nblocks = 0;
243  for (records_t::iterator it = _records.begin(); it != _records.end(); ++it)
244  {
245  record &rec = *it;
246 
247  if (Verbosity() > 0)
248  {
249  cout << "PHG4DSTReader::build_tree - Add " << rec._name << endl;
250  }
251  const string name_cnt = "n_" + rec._name;
252  const string name_cnt_desc = name_cnt + "/I";
253  _T->Branch(name_cnt.c_str(), &(rec._cnt), name_cnt_desc.c_str(),
254  BUFFER_SIZE);
255  _T->Branch(rec._name.c_str(), &(rec._arr_ptr), BUFFER_SIZE, 99);
256 
257  nblocks++;
258  }
259 
260  if (Verbosity() > 0)
261  {
262  cout << "PHG4DSTReader::build_tree - added " << nblocks << " nodes" << endl;
263  }
264  _T->SetAutoSave(16000);
265 }
266 
267 int PHG4DSTReader::process_event(PHCompositeNode *topNode)
268 {
269  // const double significand = _event / TMath::Power(10, (int) (log10(_event)));
270  //
271  // if (fmod(significand, 1.0) == 0 && significand <= 10)
272  // cout << "PHG4DSTReader::process_event - " << _event << endl;
273  _event++;
274 
275  //clean ups
276  _particle_set.clear();
277  _vertex_set.clear();
278 
279  PHG4TruthInfoContainer *truthInfoList = findNode::getClass<
280  PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
281  if (!truthInfoList)
282  {
283  if (_event < 2)
284  cout
285  << "PHG4DSTReader::process_event - Error - can not find node G4TruthInfo. Quit processing!"
286  << endl;
287  return Fun4AllReturnCodes::ABORTEVENT;
288  }
289 
290  for (records_t::iterator it = _records.begin(); it != _records.end(); ++it)
291  {
292  record &rec = *it;
293 
294  rec._cnt = 0;
295  assert(rec._arr.get() == rec._arr_ptr);
296  assert(rec._arr.get());
297  rec._arr->Clear();
298 
299  if (rec._type == record::typ_hit)
300  {
301  if (Verbosity() >= 2)
302  cout << "PHG4DSTReader::process_event - processing " << rec._name
303  << endl;
304 
305  PHG4HitContainer *hits = findNode::getClass<PHG4HitContainer>(topNode,
306  rec._name);
307  if (!hits)
308  {
309  if (_event < 2)
310  cout
311  << "PHG4DSTReader::process_event - Error - can not find node "
312  << rec._name << endl;
313  }
314  else
315  {
316  PHG4HitContainer::ConstRange hit_range = hits->getHits();
317 
318  if (Verbosity() >= 2)
319  cout << "PHG4DSTReader::process_event - processing "
320  << rec._name << " and received " << hits->size() << " hits"
321  << endl;
322 
323  for (PHG4HitContainer::ConstIterator hit_iter = hit_range.first;
324  hit_iter != hit_range.second; hit_iter++)
325  {
326  PHG4Hit *hit = hit_iter->second;
327 
328  // hit_type * hit = dynamic_cast<hit_type *>(hit_raw);
329  //
330  assert(hit);
331 
332  new ((*(rec._arr.get()))[rec._cnt]) hit_type(hit);
333 
334  hit_type *new_hit =
335  dynamic_cast<hit_type *>(rec._arr.get()->At(rec._cnt));
336  assert(new_hit);
337 
338  // for (int i = 0; i < 2; i++)
339  // {
340  // new_hit->set_x(i, hit->get_x(i));
341  // new_hit->set_y(i, hit->get_y(i));
342  // new_hit->set_z(i, hit->get_z(i));
343  // new_hit->set_t(i, hit->get_t(i));
344  // }
345  // new_hit->set_edep(hit->get_edep());
346  // new_hit->set_layer(hit->get_layer());
347  // new_hit->set_hit_id(hit->get_hit_id());
348  // new_hit->set_trkid(hit->get_trkid());
349 
350  // *new_hit = (*hit);
351 
352  if (_load_active_particle)
353  _particle_set.insert(hit->get_trkid());
354 
355  if (Verbosity() >= 2)
356  cout << "PHG4DSTReader::process_event - processing "
357  << rec._name << " and hit " << hit->get_hit_id()
358  << " with track id " << hit->get_trkid() << endl;
359 
360  rec._cnt++;
361  }
362  } // if (!hits)
363  } // if (rec._type == record::typ_hit)
364  else if (rec._type == record::typ_tower)
365  {
366  if (Verbosity() >= 2)
367  cout << "PHG4DSTReader::process_event - processing tower "
368  << rec._name << endl;
369 
370  RawTowerContainer *hits = findNode::getClass<RawTowerContainer>(
371  topNode, rec._name);
372  if (!hits)
373  {
374  if (_event < 2)
375  cout
376  << "PHG4DSTReader::process_event - Error - can not find node "
377  << rec._name << endl;
378  }
379  else
380  {
381  RawTowerContainer::ConstRange hit_range = hits->getTowers();
382 
383  if (Verbosity() >= 2)
384  cout << "PHG4DSTReader::process_event - processing "
385  << rec._name << " and received " << hits->size()
386  << " tower hits" << endl;
387 
388  for (RawTowerContainer::ConstIterator hit_iter = hit_range.first;
389  hit_iter != hit_range.second; hit_iter++)
390  {
391  RawTower *hit_raw = hit_iter->second;
392 
393  RawTower_type *hit = dynamic_cast<RawTower_type *>(hit_raw);
394  // RawTower * hit = hit_iter->second;
395 
396  assert(hit);
397 
398  if (hit->get_energy() < _tower_zero_sup)
399  {
400  if (Verbosity() >= 2)
401  cout
402  << "PHG4DSTReader::process_event - suppress low energy tower hit "
403  << rec._name << " @ ("
404  // << hit->get_thetaMin()
405  // << ", " << hit->get_phiMin()
406  << "), Energy = " << hit->get_energy() << endl;
407 
408  continue;
409  }
410 
411  new ((*(rec._arr.get()))[rec._cnt]) RawTower_type();
412 
413  if (Verbosity() >= 2)
414  cout
415  << "PHG4DSTReader::process_event - processing Tower hit "
416  << rec._name << " @ ("
417  // << hit->get_thetaMin() << ", "
418  // << hit->get_phiMin()
419  << "), Energy = " << hit->get_energy() << " - "
420  << rec._arr.get()->At(rec._cnt)->ClassName() << endl;
421 
422  // rec._arr->Print();
423 
424  RawTower_type *new_hit =
425  dynamic_cast<RawTower_type *>(rec._arr.get()->At(rec._cnt));
426  assert(new_hit);
427 
428  // //v1 copy
429  // new_hit->get_bineta(hit->get_bineta());
430  // new_hit->get_binphi(hit->get_binphi());
431  //
432  // //v2 copy
433  // new_hit->set_thetaMin(hit->get_thetaMin());
434  // new_hit->set_thetaSize(hit->get_thetaSize());
435  // new_hit->set_phiMin(hit->get_phiMin());
436  // new_hit->set_phiSize(hit->get_phiSize());
437  // new_hit->set_zMin(hit->get_zMin());
438  // new_hit->set_zSize(hit->get_zSize());
439 
440  *new_hit = (*hit);
441 
442  rec._cnt++;
443  }
444  } // if (!hits)
445  } // if (rec._type == record::typ_hit)
446  else if (rec._type == record::typ_jets)
447  {
448  // std::cout
449  // << "PHG4DSTReader::AddJet - Error - temp. disabled until jet added back to sPHENIX software"
450  // << std::endl;
451  //
452  if (Verbosity() >= 2)
453  cout << "PHG4DSTReader::process_event - processing jets "
454  << rec._name << endl;
455 
456  JetMap *hits = findNode::getClass<JetMap>(topNode, rec._name);
457  if (!hits)
458  {
459  if (_event < 2)
460  cout
461  << "PHG4DSTReader::process_event - Error - can not find node "
462  << rec._name << endl;
463  }
464  else
465  {
466  if (Verbosity() >= 2)
467  cout << "PHG4DSTReader::process_event - processing "
468  << rec._name << " and received " << hits->size() << " jets"
469  << endl;
470 
471  // for every recojet
472  for (JetMap::Iter iter = hits->begin(); iter != hits->end();
473  ++iter)
474  {
475  Jet *hit_raw = iter->second;
476 
477  if (Verbosity() >= 2)
478  cout << "PHG4DSTReader::process_event - processing jet "
479  << rec._name << " @ (" << hit_raw->get_eta() << ", "
480  << hit_raw->get_phi() << "), pT = " << hit_raw->get_pt()
481  << " - with raw type " << hit_raw->ClassName() << endl;
482 
483  PHPyJet_type *hit = dynamic_cast<PHPyJet_type *>(hit_raw);
484 
485  assert(hit);
486 
487  new ((*(rec._arr.get()))[rec._cnt]) PHPyJet_type();
488 
489  PHPyJet_type *new_hit =
490  dynamic_cast<PHPyJet_type *>(rec._arr.get()->At(rec._cnt));
491  assert(new_hit);
492 
493  *new_hit = (*hit);
494 
495  rec._cnt++;
496  }
497  } // if (!hits)
498  } // if (rec._type == record::typ_hit)
499  else if (rec._type == record::typ_part)
500  {
501  if (_load_all_particle)
502  {
503  static bool once = true;
504  if (once)
505  {
506  cout
507  << "PHG4DSTReader::process_event - will load all particle from G4TruthInfo"
508  << endl;
509 
510  once = false;
511  }
512 
513  for (auto particle_iter = truthInfoList->GetMap().begin();
514  particle_iter != truthInfoList->GetMap().end();
515  ++particle_iter)
516  {
517  _particle_set.insert(particle_iter->first);
518  }
519 
520  } // if (_load_all_particle)
521  else
522  {
523  static bool once = true;
524  if (once)
525  {
526  if (Verbosity() > 0)
527  {
528  cout << "PHG4DSTReader::process_event - will load primary particle from G4TruthInfo" << endl;
529  }
530  once = false;
531  }
532 
533  PHG4TruthInfoContainer::ConstRange primary_range =
534  truthInfoList->GetPrimaryParticleRange();
535 
536  for (PHG4TruthInfoContainer::ConstIterator particle_iter =
537  primary_range.first;
538  particle_iter != primary_range.second;
539  ++particle_iter)
540  {
541  _particle_set.insert(particle_iter->first);
542 
543  } // if (_load_all_particle) else
544  }
545  for (PartSet_t::const_iterator i = _particle_set.begin();
546  i != _particle_set.end(); ++i)
547  {
548  auto particle = truthInfoList->GetParticle(*i);
549  if (!particle)
550  {
551  cout
552  << "PHG4DSTReader::process_event - ERROR - can not find particle/track ID "
553  << *i << " in G4TruthInfo" << endl;
554 
555  continue;
556  }
557 
558  add_particle(rec, particle);
559  } // for(PartSet_t::const_iterator i = _particle_set.begin();i!=_particle_set.end();i++)
560 
561  } // if (rec._type == record::typ_part)
562  else if (rec._type == record::typ_vertex)
563  {
564  static bool once = true;
565  if (once)
566  {
567  if (Verbosity() > 0)
568  {
569  cout << "PHG4DSTReader::process_event - will load vertex from G4TruthInfo" << endl;
570  }
571  once = false;
572  }
573 
574  for (PartSet_t::const_iterator i = _vertex_set.begin();
575  i != _vertex_set.end(); ++i)
576  {
577  PHG4VtxPoint *v = truthInfoList->GetVtx(*i);
578  if (!v)
579  {
580  cout
581  << "PHG4DSTReader::process_event - ERROR - can not find vertex ID "
582  << *i << " in G4TruthInfo" << endl;
583 
584  continue;
585  }
586 
587  new ((*(rec._arr.get()))[rec._cnt]) vertex_type();
588 
589  if (Verbosity() >= 2)
590  cout << "PHG4DSTReader::process_event - saving vertex id "
591  << *i << endl;
592 
593  vertex_type *new_v =
594  static_cast<vertex_type *>(rec._arr.get()->At(rec._cnt));
595  assert(new_v);
596 
597  new_v->set_x(v->get_x());
598  new_v->set_y(v->get_y());
599  new_v->set_z(v->get_z());
600  new_v->set_t(v->get_t());
601  new_v->set_id(v->get_id());
602 
603  rec._cnt++;
604  } // else if (rec._type == record::typ_vertex)
605 
606  } // if (_load_all_particle)
607 
608  } // for (records_t::iterator it = _records.begin(); it != _records.end(); ++it)
609 
610  if (_T)
611  _T->Fill();
612 
613  return 0;
614 } // for (records_t::iterator it = _records.begin(); it != _records.end(); ++it)
615 
616 void PHG4DSTReader::add_particle(PHG4DSTReader::record &rec, PHG4Particle *part)
617 {
618  assert(part);
619 
620  // if (Verbosity() >= 2)
621  // cout << "PHG4DSTReader::process_event - Particle type is "
622  // << p_raw->GetName() << " - " << p_raw->ClassName()
623  // << " get_track_id = " << p_raw->get_track_id() << endl;
624 
625  // part_type * part = dynamic_cast<part_type *>(p_raw);
626 
627  // assert(part);
628 
629  new ((*(rec._arr.get()))[rec._cnt]) part_type();
630 
631  part_type *new_part = static_cast<part_type *>(rec._arr.get()->At(rec._cnt));
632  assert(new_part);
633 
634  new_part->set_track_id(part->get_track_id());
635  new_part->set_vtx_id(part->get_vtx_id());
636  new_part->set_parent_id(part->get_parent_id());
637  new_part->set_primary_id(part->get_primary_id());
638  new_part->set_name(part->get_name());
639  new_part->set_pid(part->get_pid());
640  new_part->set_px(part->get_px());
641  new_part->set_py(part->get_py());
642  new_part->set_pz(part->get_pz());
643  new_part->set_e(part->get_e());
644 
645  _vertex_set.insert(part->get_vtx_id());
646 
647  rec._cnt++;
648 }
649 
650 int PHG4DSTReader::End(PHCompositeNode * /*topNode*/)
651 {
652  if (Verbosity() > 1)
653  {
654  cout << "PHG4DSTReader::End - Clean ups" << endl;
655  }
656 
657  if (_T)
658  {
659  PHTFileServer::get().cd(_out_file_name);
660  _T->Write();
661  _T->ResetBranchAddresses();
662  }
663 
664  _records.clear();
665 
666  return 0;
667 }