ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4TowerNtuple.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4TowerNtuple.cc
1 #include "G4TowerNtuple.h"
2 
3 #include <calobase/RawTower.h>
4 #include <calobase/RawTowerContainer.h>
5 #include <calobase/RawTowerGeomContainer.h>
6 
8 #include <fun4all/SubsysReco.h> // for SubsysReco
9 
10 #include <phool/getClass.h>
11 
12 #include <TFile.h>
13 #include <TH1.h>
14 #include <TNtuple.h>
15 
16 #include <cmath> // for isfinite
17 #include <iostream> // for operator<<, basic_ostream
18 #include <sstream>
19 #include <utility> // for pair, make_pair
20 
21 using namespace std;
22 
23 G4TowerNtuple::G4TowerNtuple(const std::string &name, const std::string &filename)
24  : SubsysReco(name)
25  , nblocks(0)
26  , hm(nullptr)
27  , _filename(filename)
28  , ntup(nullptr)
29  , outfile(nullptr)
30 {
31 }
32 
34 {
35  // delete ntup;
36  delete hm;
37 }
38 
40 {
41  hm = new Fun4AllHistoManager(Name());
42  outfile = new TFile(_filename.c_str(), "RECREATE");
43  ntup = new TNtuple("towerntup", "G4Towers", "detid:phi:eta:energy");
44  // ntup->SetDirectory(0);
45  TH1 *h1 = new TH1F("energy1GeV", "energy 0-1GeV", 1000, 0, 1);
46  eloss.push_back(h1);
47  h1 = new TH1F("energy100GeV", "energy 0-100GeV", 1000, 0, 100);
48  eloss.push_back(h1);
49  return 0;
50 }
51 
53 {
54  ostringstream nodename;
55  ostringstream geonodename;
56  set<string>::const_iterator iter;
57  vector<TH1 *>::const_iterator eiter;
58  for (iter = _node_postfix.begin(); iter != _node_postfix.end(); ++iter)
59  {
60  int detid = (_detid.find(*iter))->second;
61  nodename.str("");
62  nodename << "TOWER_" << _tower_type[*iter];
63  geonodename.str("");
64  geonodename << "TOWERGEOM_" << *iter;
65  RawTowerGeomContainer *towergeom = findNode::getClass<RawTowerGeomContainer>(topNode, geonodename.str());
66  if (!towergeom)
67  {
68  cout << "no geometry node " << geonodename.str() << " for " << *iter << endl;
69  continue;
70  }
71  RawTowerContainer *towers = findNode::getClass<RawTowerContainer>(topNode, nodename.str());
72  if (towers)
73  {
74  double esum = 0;
75  // double numtowers = towers->size();
76  // ntowers[i]->Fill(numtowers);
77  // cout << "number of towers: " << towers->size() << endl;
78  RawTowerContainer::ConstRange tower_range = towers->getTowers();
79  for (RawTowerContainer::ConstIterator tower_iter = tower_range.first; tower_iter != tower_range.second; tower_iter++)
80 
81  {
82  double energy = tower_iter->second->get_energy();
83  if (!isfinite(energy))
84  {
85  cout << "invalid energy: " << energy << endl;
86  }
87  esum += tower_iter->second->get_energy();
88  int phibin = tower_iter->second->get_binphi();
89  int etabin = tower_iter->second->get_bineta();
90  // to search the map fewer times, cache the geom object until the layer changes
91  double phi = towergeom->get_phicenter(phibin);
92  double eta = towergeom->get_etacenter(etabin);
93  ntup->Fill(detid,
94  phi,
95  eta,
96  tower_iter->second->get_energy());
97  }
98  for (eiter = eloss.begin(); eiter != eloss.end(); ++eiter)
99  {
100  (*eiter)->Fill(esum);
101  }
102  }
103  }
104  return 0;
105 }
106 
108 {
109  outfile->cd();
110  ntup->Write();
111  outfile->Write();
112  outfile->Close();
113  delete outfile;
114  hm->dumpHistos(_filename, "UPDATE");
115  return 0;
116 }
117 
118 void G4TowerNtuple::AddNode(const std::string &name, const std::string &twrtype, const int detid)
119 {
120  ostringstream twrname;
121  twrname << twrtype << "_" << name;
122  _node_postfix.insert(name);
123  _tower_type.insert(make_pair(name, twrname.str()));
124  _detid[name] = detid;
125  return;
126 }