Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Calib.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file Calib.cc
1 #include "Calib.h"
2 #include <fun4all/Fun4AllReturnCodes.h>
3 #include <fun4all/PHTFileServer.h>
4 #include <fun4all/Fun4AllServer.h>
5 
6 #include <g4main/PHG4HitContainer.h>
7 #include <g4main/PHG4Hit.h>
8 #include <g4main/PHG4TruthInfoContainer.h>
9 #include <g4main/PHG4VtxPoint.h>
10 #include <g4main/PHG4Particle.h>
11 
12 #include <phool/PHCompositeNode.h>
13 #include <phool/PHIODataNode.h>
14 #include <phool/PHNode.h>
15 #include <phool/PHNodeIterator.h>
16 #include <phool/PHObject.h>
17 #include <phool/getClass.h>
18 #include <phool/phool.h>
19 
20 
21 #include <calobase/RawTowerDefs.h>
22 #include <calobase/RawTowerContainer.h>
23 #include <calobase/RawTower.h>
24 #include <calobase/RawTowerGeomContainer.h>
25 #include <calobase/RawTowerGeom.h>
26 #include <calobase/RawClusterContainer.h>
27 #include <calobase/RawCluster.h>
28 #include <calobase/RawTowerv1.h>
29 
30 #include <phgeom/PHGeomUtility.h>
31 #include <TTree.h>
32 #include <TH2D.h>
33 #include <TProfile2D.h>
34 #include <TVector3.h>
35 #include <TRandom3.h>
36 #include <TMath.h>
37 
38 #include <fstream>
39 
40 
41 
42 using namespace std;
43 
44 //----------------------------------------------------------------------------//
45 //-- Constructor:
46 //-- simple initialization
47 //----------------------------------------------------------------------------//
48 
49 Calib::Calib(const string &name) :
50  SubsysReco(name), _event_tree(NULL)
51 {
52  //initialize
53  _event = 0;
54  _outfile_name = "Tower_test.root";
55 
56 }
57 
58 //----------------------------------------------------------------------------//
59 //-- Init():
60 //-- Intialize all histograms, trees, and ntuples
61 //----------------------------------------------------------------------------//
62 int Calib::Init(PHCompositeNode *topNode) {
63 
64  cout << PHWHERE << " Opening file " << _outfile_name << endl;
65 
66  PHTFileServer::get().open(_outfile_name, "RECREATE");
67  PHTFileServer::get().cd(_outfile_name);
68 
69  Double_t xbinlimits[25] =
70  {
71  -1.10001,-1.00833,-0.916665,-0.824998,-0.733332,-0.641665,-0.549998,-0.458332,-0.366665,-0.274998,-0.183332,-0.091665,0.,0.0916683,0.183335,0.275002,0.366668,0.458335,0.550002,0.641668,0.733335,0.825002,0.916668,1.00833,1.10001
72  };
73 
74  Double_t ybinlimits[65] =
75  {
76  -3.16565,-3.01936,-2.92119,-2.82301,-2.72484,-2.62666,-2.52849,-2.43032,-2.33214,-2.23397,-2.13579,-2.03762,-1.93944,-1.84127,-1.74309,-1.64492,-1.54674,-1.44857,-1.35039,-1.25222,-1.15404,-1.05587,-0.95769,-0.859519,-0.761344,-0.663169,-0.564994,-0.46682,-0.368645,-0.27047,-0.172295,-0.074121,0.024068900,0.122229,0.220404,0.318578,0.416753,0.514928,0.613103,0.711278,0.809452,0.907627,1.0058,1.10398,1.20215,1.30033,1.3985,1.49668,1.59485,1.69303,1.7912,1.88937,1.98755,2.08572,2.1839,2.28207,2.38025,2.47842,2.5766,2.67477,2.77295,2.87112,2.9693,3.06747,3.16565
77  };
78 
79  Long_t xbinnumb = sizeof(xbinlimits)/sizeof(Double_t) - 1;
80  Long_t ybinnumb = sizeof(ybinlimits)/sizeof(Double_t) - 1;
81 
82 
83  hprof2d = new TProfile2D("hprof2d","Profile of tower energy versus eta and phi", xbinnumb,xbinlimits, ybinnumb,ybinlimits);
84  hprof2d->SetXTitle("#eta");
85  hprof2d->SetYTitle("#phi");
86 
87  _event_tree = new TTree("event", "Calib => event info");
88 
89  _event_tree->Branch("event", &_event, "_event/I");
90  _event_tree->Branch("oCalib_ohcal", &_oCalib_ohcal, "_oCalib_ohcal/F");
91 
92 
93  return Fun4AllReturnCodes::EVENT_OK;
94 }
95 
96 int Calib::InitRun(PHCompositeNode *topNode)
97 {
98  return Fun4AllReturnCodes::EVENT_OK;
99 }
100 
101 //----------------------------------------------------------------------------//
102 //-- process_event():
103 //-- Call user instructions for every event.
104 //-- This function contains the analysis structure.
105 //----------------------------------------------------------------------------//
106 
107 int Calib::process_event(PHCompositeNode *topNode)
108 {
109  _event++;
110 
111  GetNodes(topNode);
112 
113  fill_tree(topNode);
114 
115  return Fun4AllReturnCodes::EVENT_OK;
116 }
117 
118 //----------------------------------------------------------------------------//
119 //-- End():
120 //-- End method, wrap everything up
121 //----------------------------------------------------------------------------//
122 
123 int Calib::EndRun(PHCompositeNode *topNode)
124 {
125 
126  return Fun4AllReturnCodes::EVENT_OK;
127 
128 }
129 
130 int Calib::End(PHCompositeNode *topNode)
131 {
132 
133  PHTFileServer::get().cd(_outfile_name);
134  PHTFileServer::get().write(_outfile_name);
135 
136  return Fun4AllReturnCodes::EVENT_OK;
137 }
138 
139 //----------------------------------------------------------------------------//
140 //-- fill_tree():
141 //-- Fill the various trees...
142 //----------------------------------------------------------------------------//
143 
144 void Calib::fill_tree(PHCompositeNode *topNode)
145 {
146 
147  cout << _event << endl;
148 
149 
150 
151  //ohcal towers
152  int ietabin, iphibin = -1;
153  float oh_ocal = 0.0;
154  float OHCAL_o[24][64] = {{0.0}};
155  float OHCAL_tower_eta[24] = {0.0};
156  float OHCAL_tower_phi[64] = {0.0};
157 
158  //loop over ohcal towers
159  RawTowerContainer::ConstRange begin_end_oo = _ohcal_towers_o->getTowers();
160  RawTowerContainer::ConstIterator oh_o = begin_end_oo.first;
161  for (; oh_o != begin_end_oo.second; ++oh_o) {
162  RawTowerDefs::keytype towerid_o_o = oh_o->first;
163  RawTower *rawtower_o_o = _ohcal_towers_o->getTower(towerid_o_o);
164  if(rawtower_o_o) {
165  RawTowerGeom *tgeo_o_o = _ohcal_towergeom->get_tower_geometry(towerid_o_o);
166  // get towers energy
167  oh_ocal = rawtower_o_o->get_energy();
168 
169  // binning in eta and phi
170  ietabin = tgeo_o_o->get_column();
171  iphibin = tgeo_o_o->get_row();
172  if((iphibin >= 0)&&(ietabin >= 0))
173  {
174  OHCAL_o[ietabin][iphibin] += oh_ocal;
175  OHCAL_tower_eta[ietabin] = tgeo_o_o->get_eta();
176  OHCAL_tower_phi[iphibin] = tgeo_o_o->get_phi();
177  }
178 
179  }
180  }
181 
182 
183 
184  //---- Store energy ------//
185  _oCalib_ohcal = 0.0;
186 
187  for(int i=0; i<24; i++){
188  for(int j=0; j<64; j++){
189  _oCalib_ohcal += OHCAL_o[i][j];
190  hprof2d->Fill(OHCAL_tower_eta[i], OHCAL_tower_phi[j],OHCAL_o[i][j]);
191  }
192  }
193 
194 
195 _event_tree->Fill();
196 
197 
198 return;
199 
200 }
201 
202 //----------------------------------------------------------------------------//
203 //-- GetNodes():
204 //-- Get all the all the required nodes off the node tree
205 //----------------------------------------------------------------------------//
206 int Calib::GetNodes(PHCompositeNode * topNode) {
207 
208  //oHCAL
209 
210  _ohcal_towers_o = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_HCALOUT");
211  if (!_ohcal_towers_o)
212  {
213  std::cout << PHWHERE << ": Could not find node TOWER_CALIB_HCALOUT" << std::endl;
214  return Fun4AllReturnCodes::ABORTEVENT;
215  }
216 
217 
218  //towergeom
219  _ohcal_towergeom = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALOUT");
220  if (! _ohcal_towergeom)
221  {
222  cout << PHWHERE << ": Could not find node TOWERGEOM_HCALOUT" << endl;
223  return Fun4AllReturnCodes::ABORTEVENT;
224  }
225 
226  return Fun4AllReturnCodes::EVENT_OK;
227 }
228