Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
HCalJetPhiShift.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file HCalJetPhiShift.cc
1 //_________________________________________________________________________________..
2 // Veronica Verkest, June 2023
3 // module for performing an analysis on the phi shift between truth and calo jets
4 //_________________________________________________________________________________..
5 
6 #include "HCalJetPhiShift.h"
7 
8 #include "g4eval/CaloEvalStack.h"
9 #include "g4eval/CaloRawClusterEval.h"
10 #include "g4eval/CaloRawTowerEval.h"
11 #include "g4eval/CaloTruthEval.h"
12 
13 #include <g4main/PHG4Particle.h>
14 #include <g4main/PHG4Shower.h>
15 #include <g4main/PHG4TruthInfoContainer.h>
16 #include <g4main/PHG4VtxPoint.h>
17 
18 //#include <g4detectors/PHG4CellContainer.h>
19 //#include <g4detectors/PHG4CylinderCellGeomContainer.h>
20 //#include <g4detectors/PHG4CylinderGeomContainer.h>
21 //#include <g4detectors/PHG4Cell.h>
22 //#include <g4detectors/PHG4CylinderCellGeom.h>
23 
24 #include <fun4all/Fun4AllReturnCodes.h>
25 #include <fun4all/PHTFileServer.h>
26 
27 #include <phool/PHCompositeNode.h>
28 #include <phool/getClass.h>
29 
30 #include <calobase/RawTower.h>
31 #include <calobase/RawTowerContainer.h>
32 #include <calobase/RawTowerGeom.h>
33 #include <calobase/RawTowerGeomContainer.h>
34 #include <calobase/TowerInfoContainer.h>
35 #include <calobase/TowerInfo.h>
36 
37 #include <TFile.h>
38 #include <TTree.h>
39 
40 //____________________________________________________________________________..
41 HCalJetPhiShift::HCalJetPhiShift(const std::string &name, const std::string &outputFile):
42 SubsysReco(name),
43 m_outputFileName(outputFile),
44 m_T(nullptr),
45 m_event(-1),
46 m_nTow_in(0),
47 m_nTow_out(0),
48 m_nTow_emc(0),
49 m_eta(),
50 m_phi(),
51 m_e(),
52 m_pt(),
53 m_vx(),
54 m_vy(),
55 m_vz(),
56 m_id(),
57 m_eta_in(),
58 m_phi_in(),
59 m_e_in(),
60 m_ieta_in(),
61 m_iphi_in(),
62 m_eta_out(),
63 m_phi_out(),
64 m_e_out(),
65 m_ieta_out(),
66 m_iphi_out(),
67 m_eta_emc(),
68 m_phi_emc(),
69 m_e_emc(),
70 m_ieta_emc(),
71 m_iphi_emc()
72 
73 {
74  std::cout << "HCalJetPhiShift::HCalJetPhiShift(const std::string &name) Calling ctor" << std::endl;
75 }
76 
77 //____________________________________________________________________________..
78 HCalJetPhiShift::~HCalJetPhiShift()
79 {
80  std::cout << "HCalJetPhiShift::~HCalJetPhiShift() Calling dtor" << std::endl;
81 }
82 
83 //____________________________________________________________________________..
84 int HCalJetPhiShift::Init(PHCompositeNode* /*topNode*/)
85 {
86  std::cout << "HCalJetPhiShift::Init(PHCompositeNode *topNode) Initializing" << std::endl;
87  PHTFileServer::get().open(m_outputFileName, "RECREATE");
88  std::cout << "HCalJetPhiShift::Init - Output to " << m_outputFileName << std::endl;
89 
90  // configure Tree
91  m_T = new TTree("T", "HCalJetPhiShift Tree");
92  m_T->Branch("event", &m_event);
93  m_T->Branch("eta", &m_eta);
94  m_T->Branch("phi", &m_phi);
95  m_T->Branch("e", &m_e);
96  m_T->Branch("pt", &m_pt);
97  m_T->Branch("vx", &m_vx);
98  m_T->Branch("vy", &m_vy);
99  m_T->Branch("vz", &m_vz);
100  m_T->Branch("id", &m_id);
101  m_T->Branch("nTow_in", &m_nTow_in);
102  m_T->Branch("eta_in", &m_eta_in);
103  m_T->Branch("phi_in", &m_phi_in);
104  m_T->Branch("e_in", &m_e_in);
105  m_T->Branch("ieta_in", &m_ieta_in);
106  m_T->Branch("iphi_in", &m_iphi_in);
107  m_T->Branch("nTow_out", &m_nTow_out);
108  m_T->Branch("eta_out", &m_eta_out);
109  m_T->Branch("phi_out", &m_phi_out);
110  m_T->Branch("e_out", &m_e_out);
111  m_T->Branch("ieta_out", &m_ieta_out);
112  m_T->Branch("iphi_out", &m_iphi_out);
113  m_T->Branch("nTow_emc", &m_nTow_emc);
114  m_T->Branch("eta_emc", &m_eta_emc);
115  m_T->Branch("phi_emc", &m_phi_emc);
116  m_T->Branch("e_emc", &m_e_emc);
117  m_T->Branch("ieta_emc", &m_ieta_emc);
118  m_T->Branch("iphi_emc", &m_iphi_emc);
119 
120 
121  return Fun4AllReturnCodes::EVENT_OK;
122 }
123 
124 //____________________________________________________________________________..
125 int HCalJetPhiShift::InitRun(PHCompositeNode *topNode)
126 {
127  std::cout << "HCalJetPhiShift::InitRun(PHCompositeNode *topNode) Initializing for Run XXX" << std::endl;
128  return Fun4AllReturnCodes::EVENT_OK;
129 }
130 
131 //____________________________________________________________________________..
132 int HCalJetPhiShift::process_event(PHCompositeNode *topNode)
133 {
134  // std::cout << "HCalJetPhiShift::process_event(PHCompositeNode *topNode) Processing Event" << std::endl;
135  ++m_event;
136 
137  //fill the tree
138  FillTTree(topNode);
139 
140  return Fun4AllReturnCodes::EVENT_OK;
141 }
142 
143 //____________________________________________________________________________..
144 int HCalJetPhiShift::ResetEvent(PHCompositeNode *topNode)
145 {
146  // std::cout << "HCalJetPhiShift::ResetEvent(PHCompositeNode *topNode) Resetting internal structures, prepare for next event" << std::endl;
147  m_nTow_in = 0;
148  m_nTow_out = 0;
149  m_nTow_emc = 0;
150  m_id.clear();
151  m_eta_in.clear();
152  m_phi_in.clear();
153  m_e_in.clear();
154  m_ieta_in.clear();
155  m_iphi_in.clear();
156 
157  m_eta_out.clear();
158  m_phi_out.clear();
159  m_e_out.clear();
160  m_ieta_out.clear();
161  m_iphi_out.clear();
162 
163  m_eta_emc.clear();
164  m_phi_emc.clear();
165  m_e_emc.clear();
166  m_ieta_emc.clear();
167  m_iphi_emc.clear();
168 
169  return Fun4AllReturnCodes::EVENT_OK;
170 }
171 
172 //____________________________________________________________________________..
173 int HCalJetPhiShift::EndRun(const int runnumber)
174 {
175  std::cout << "HCalJetPhiShift::EndRun(const int runnumber) Ending Run for Run " << runnumber << std::endl;
176  return Fun4AllReturnCodes::EVENT_OK;
177 }
178 
179 //____________________________________________________________________________..
180 int HCalJetPhiShift::End(PHCompositeNode *topNode)
181 {
182  std::cout << "HCalJetPhiShift::End - Output to " << m_outputFileName << std::endl;
183  PHTFileServer::get().cd(m_outputFileName);
184 
185  m_T->Write();
186  std::cout << "HCalJetPhiShift::End(PHCompositeNode *topNode) This is the End..." << std::endl;
187  return Fun4AllReturnCodes::EVENT_OK;
188 }
189 
190 //____________________________________________________________________________..
191 int HCalJetPhiShift::Reset(PHCompositeNode *topNode)
192 {
193  std::cout << "HCalJetPhiShift::Reset(PHCompositeNode *topNode) being Reset" << std::endl;
194  return Fun4AllReturnCodes::EVENT_OK;
195 }
196 
197 //____________________________________________________________________________..
198 void HCalJetPhiShift::Print(const std::string &what) const
199 {
200  std::cout << "HCalJetPhiShift::Print(const std::string &what) const Printing info for " << what << std::endl;
201 }
202 
203 //____________________________________________________________________________..
204 int HCalJetPhiShift::FillTTree(PHCompositeNode *topNode)
205 {
206 // std::cout << "HCalJetPhiShift::FillTTree(PHCompositeNode *topNode)" << std::endl;
207 
208  // fill truth info from nodetree
209  PHG4TruthInfoContainer* truthinfo = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
210  if (!truthinfo)
211  {
212  std::cout << PHWHERE << " ERROR: Can't find G4TruthInfo" << std::endl;
213  exit(-1);
214  }
215 
216  PHG4TruthInfoContainer::ConstRange range = truthinfo->GetPrimaryParticleRange();
217  for (PHG4TruthInfoContainer::ConstIterator iter = range.first;
218  iter != range.second;
219  ++iter)
220  {
221  PHG4Particle* primary = iter->second;
222 
223  if (primary->get_e() < 0.)
224  {
225  continue;
226  }
227 
228  float gpx = primary->get_px();
229  float gpy = primary->get_py();
230  float gpz = primary->get_pz();
231  m_e = primary->get_e();
232 
233  m_pt = std::sqrt(gpx * gpx + gpy * gpy);
234  m_eta = NAN;
235  if (m_pt != 0.0)
236  {
237  m_eta = std::asinh(gpz / m_pt);
238  }
239  m_phi = std::atan2(gpy, gpx);
240 
241  PHG4VtxPoint* gvertex = truthinfo->GetPrimaryVtx(truthinfo->GetPrimaryVertexIndex());
242  m_vx = gvertex->get_x();
243  m_vy = gvertex->get_y();
244  m_vz = gvertex->get_z();
245  }
246 
247  //calorimeter towers
248  TowerInfoContainer *towersIH3 = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_HCALIN");
249  TowerInfoContainer *towersOH3 = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_HCALOUT");
250 // TowerInfoContainer *towersEM3 = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_CEMC");
251  TowerInfoContainer *towersEM3 = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_RAW_CEMC");
252  RawTowerGeomContainer *tower_geomIH = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALIN");
253  RawTowerGeomContainer *tower_geomOH = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALOUT");
254  RawTowerGeomContainer *tower_geomEM = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_CEMC");
255  if(!towersIH3 || !towersOH3 || !towersEM3){
256  std::cout
257  <<"HCalJetPhiShift::process_event - Error cannot find raw tower node "
258 
259  << std::endl;
260  exit(-1);
261  }
262 
263  if(!tower_geomIH || !tower_geomOH || !tower_geomEM){
264  std::cout
265  <<"HCalJetPhiShift::process_event - Error cannot find raw tower geometry "
266 
267  << std::endl;
268  exit(-1);
269  }
270 
271  TowerInfo *tower_in, *tower_out;//, *tower_emc;
272 
273  const int n_channels_IH = (int) towersIH3->size();
274 
275  // Inner HCal
276  for (int i_chan=0; i_chan<n_channels_IH; ++i_chan) {
277  tower_in = towersIH3->get_tower_at_channel(i_chan);
278  if (tower_in->get_energy()>0.) {
279  ++m_nTow_in;
280 
281  unsigned int calokey = towersIH3->encode_key(i_chan);
282  int ieta = towersIH3->getTowerEtaBin(calokey);
283  int iphi = towersIH3->getTowerPhiBin(calokey);
284 
285  const RawTowerDefs::keytype key = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::HCALIN, ieta, iphi);
286  float tower_phi = tower_geomIH->get_tower_geometry(key)->get_phi();
287  float tower_eta = tower_geomIH->get_tower_geometry(key)->get_eta();
288 
289  m_eta_in.push_back(tower_eta);
290  m_phi_in.push_back(tower_phi);
291  m_e_in.push_back(tower_in->get_energy());
292  m_ieta_in.push_back(ieta);
293  m_iphi_in.push_back(iphi);
294  }
295 
296  // Outer HCal
297  tower_out = towersOH3->get_tower_at_channel(i_chan);
298  if (tower_out->get_energy()>0.) {
299  ++m_nTow_out;
300 
301  unsigned int calokey = towersOH3->encode_key(i_chan);
302  int ieta = towersOH3->getTowerEtaBin(calokey);
303  int iphi = towersOH3->getTowerPhiBin(calokey);
304 
305  const RawTowerDefs::keytype key = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::HCALOUT, ieta, iphi);
306  float tower_phi = tower_geomOH->get_tower_geometry(key)->get_phi();
307  float tower_eta = tower_geomOH->get_tower_geometry(key)->get_eta();
308 
309  m_eta_out.push_back(tower_eta);
310  m_phi_out.push_back(tower_phi);
311  m_e_out.push_back(tower_out->get_energy());
312  m_ieta_out.push_back(ieta);
313  m_iphi_out.push_back(iphi);
314  }
315 
316  }
317 
318  const int n_channels_EM = (int) towersEM3->size();
319  for (int i_chan=0; i_chan<n_channels_EM; ++i_chan) {
320  // EMCal
321 // tower_emc = towersEM3->get_tower_at_channel(i_chan);
322  TowerInfo *tower_emc = towersEM3->get_tower_at_channel(i_chan);
323  if (tower_emc->get_energy()>0.) {
324  ++m_nTow_emc;
325 
326  unsigned int calokey = towersEM3->encode_key(i_chan);
327  int ieta = towersEM3->getTowerEtaBin(calokey);
328  int iphi = towersEM3->getTowerPhiBin(calokey);
329 
330  const RawTowerDefs::keytype key = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::CEMC, ieta, iphi);
331  RawTowerGeom *tower_geom = tower_geomEM->get_tower_geometry(key);
332  float tower_phi = tower_geom->get_phi();
333  float tower_eta = tower_geom->get_eta();
334 
335 // unsigned int calokey = towersEM3->encode_key(i_chan);
336 // int ieta = towersEM3->getTowerEtaBin(calokey);
337 // int iphi = towersEM3->getTowerPhiBin(calokey);
338 //
339 // const RawTowerDefs::keytype key = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::CEMC, ieta, iphi);
340 // float tower_phi = tower_geomEM->get_tower_geometry(key)->get_phi();
341 // float tower_eta = tower_geomEM->get_tower_geometry(key)->get_eta();
342 
343  m_eta_emc.push_back(tower_eta);
344  m_phi_emc.push_back(tower_phi);
345  m_e_emc.push_back(tower_emc->get_energy());
346  m_ieta_emc.push_back(ieta);
347  m_iphi_emc.push_back(iphi);
348  }
349 
350  }
351 
352 
353  m_T->Fill();
354 
355  return Fun4AllReturnCodes::EVENT_OK;
356 }