Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
HCALAnalysis.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file HCALAnalysis.C
1 // General sPHENIX tools
2 
3 #include "HCALAnalysis.h"
4 
5 #include <calobase/RawTowerContainer.h>
6 #include <pdbcalbase/PdbParameterMap.h>
7 #include <phparameter/PHParameters.h>
8 #include <ffaobjects/EventHeader.h>
9 #include <g4detectors/PHG4ScintillatorSlat.h>
10 #include <g4detectors/PHG4ScintillatorSlatContainer.h>
11 #include <g4detectors/PHG4ScintillatorSlatDefs.h>
12 
13 #include <fun4all/SubsysReco.h>
14 #include <fun4all/Fun4AllServer.h>
15 #include <fun4all/PHTFileServer.h>
16 #include <phool/PHCompositeNode.h>
17 #include <fun4all/Fun4AllReturnCodes.h>
18 #include <phool/getClass.h>
19 #include <fun4all/Fun4AllHistoManager.h>
20 #include <phool/getClass.h>
21 
22 
23 #include <g4main/PHG4HitContainer.h>
24 #include <g4main/PHG4Hit.h>
25 
26 // Root classes
27 #include <TH1.h>
28 #include <TH2.h>
29 #include <TPythia6.h>
30 #include <TMath.h>
31 #include <TFile.h>
32 #include <TNtuple.h>
33 #include <TLorentzVector.h>
34 
35 #include <assert.h>
36 
37 
38 HCALAnalysis::HCALAnalysis(const char *outfile) : SubsysReco("HCAL Analysis Development Code")
39 {
40  OutFileName = outfile;
41 
42  return;
43 }
44 
45 int HCALAnalysis::Init(PHCompositeNode *topNode)
46 {
47 
48  outputfile = new TFile(OutFileName,"RECREATE");
49 
50  calenergy = new TNtuple("calenergy","calenergy","e_hcin:e_hcout:e_cemc:ea_hcin:ea_hcout:ea_cemc:ev_hcin:ev_hcout:ev_cemc:e_magnet:e_bh:e_emcelect:e_hcalin_spt:sfemc:sfihcal:sfohcal:sfvemc:sfvihcal:sfvohcal:LCG:emc_LCG:hcalin_LCG:hcalout_LCG:e_svtx:e_svtx_support:e_bh_plus:e_bh_minus");
51 
52  return 0;
53 }
54 
55 int HCALAnalysis::process_event(PHCompositeNode *topNode)
56 {
57  GetNodes(topNode);
58 
59  /*
60  std::cout << " This event has " << _hcalout_hit_container->size() << " HCALOUT G4 hits" << std::endl;
61  std::cout << " This event has " << _hcalin_hit_container->size() << " HCALIN G4 hits" << std::endl;
62  std::cout << " This event has " << _cemc_hit_container->size() << " CEMC G4 hits" << std::endl;
63 
64  std::cout << " This event has " << _hcalout_abs_hit_container->size() << " ABSORBER_HCALOUT G4 hits" << std::endl;
65  std::cout << " This event has " << _hcalin_abs_hit_container->size() << " ABSORBER_HCALIN G4 hits" << std::endl;
66  std::cout << " This event has " << _cemc_abs_hit_container->size() << " ABSORBER_CEMC G4 hits" << std::endl;
67  */
68 
69  // Sum up energy in scintillator, absorber, and light yield for
70  // the three calorimeter sections.
71 
72  float e_hcin = 0.0, e_hcout = 0.0, e_cemc = 0.0;
73  float ea_hcin = 0.0, ea_hcout = 0.0, ea_cemc = 0.0;
74  float ev_hcin = 0.0, ev_hcout = 0.0, ev_cemc = 0.0;
75  float e_magnet = 0.0, e_bh = 0.0;
76  float e_emcelect = 0.0, e_hcalin_spt = 0.0;
77  float e_svtx = 0.0, e_svtx_support = 0.0;
78  float e_bh_plus = 0.0, e_bh_minus = 0.0;
79 
80  float lcg = 0.0; //longitudinal center of gravity
81  float emc_lcg = 0.0;
82  float hcalin_lcg = 0.0;
83  float hcalout_lcg = 0.0;
84 
85  PHG4HitContainer::ConstRange hcalout_hit_range = _hcalout_hit_container->getHits();
86  for (PHG4HitContainer::ConstIterator hit_iter = hcalout_hit_range.first;
87  hit_iter != hcalout_hit_range.second; hit_iter++){
88 
89  PHG4Hit *this_hit = hit_iter->second;
90  assert(this_hit);
91 
92  e_hcout += this_hit->get_edep();
93  ev_hcout += this_hit->get_light_yield();
94 
95  float rin = sqrt(pow(this_hit->get_x(0),2)+pow(this_hit->get_y(0),2));
96  float rout = sqrt(pow(this_hit->get_x(1),2)+pow(this_hit->get_y(1),2));
97 
98  lcg += this_hit->get_edep()*(rin + ((rout-rin)/2.0));
99  hcalout_lcg += this_hit->get_edep()*(rin + ((rout-rin)/2.0));
100 
101 
102  }
103 
104  PHG4HitContainer::ConstRange hcalin_hit_range = _hcalin_hit_container->getHits();
105  for (PHG4HitContainer::ConstIterator hit_iter = hcalin_hit_range.first;
106  hit_iter != hcalin_hit_range.second; hit_iter++){
107 
108  PHG4Hit *this_hit = hit_iter->second;
109  assert(this_hit);
110 
111  e_hcin += this_hit->get_edep();
112  ev_hcin += this_hit->get_light_yield();
113 
114  float rin = sqrt(pow(this_hit->get_x(0),2)+pow(this_hit->get_y(0),2));
115  float rout = sqrt(pow(this_hit->get_x(1),2)+pow(this_hit->get_y(1),2));
116 
117  lcg += this_hit->get_edep()*(rin + ((rout-rin)/2.0));
118  hcalin_lcg += this_hit->get_edep()*(rin + ((rout-rin)/2.0));
119 
120  }
121 
122  PHG4HitContainer::ConstRange cemc_hit_range = _cemc_hit_container->getHits();
123  for (PHG4HitContainer::ConstIterator hit_iter = cemc_hit_range.first;
124  hit_iter != cemc_hit_range.second; hit_iter++){
125 
126  PHG4Hit *this_hit = hit_iter->second;
127  assert(this_hit);
128 
129  e_cemc += this_hit->get_edep();
130  ev_cemc += this_hit->get_light_yield();
131 
132  float rin = sqrt(pow(this_hit->get_x(0),2)+pow(this_hit->get_y(0),2));
133  float rout = sqrt(pow(this_hit->get_x(1),2)+pow(this_hit->get_y(1),2));
134 
135  lcg += this_hit->get_edep()*(rin + ((rout-rin)/2.0));
136  emc_lcg += this_hit->get_edep()*(rin + ((rout-rin)/2.0));
137 
138  }
139 
140  PHG4HitContainer::ConstRange hcalout_abs_hit_range = _hcalout_abs_hit_container->getHits();
141  for (PHG4HitContainer::ConstIterator hit_iter = hcalout_abs_hit_range.first;
142  hit_iter != hcalout_abs_hit_range.second; hit_iter++){
143 
144  PHG4Hit *this_hit = hit_iter->second;
145  assert(this_hit);
146 
147  ea_hcout += this_hit->get_edep();
148 
149  float rin = sqrt(pow(this_hit->get_x(0),2)+pow(this_hit->get_y(0),2));
150  float rout = sqrt(pow(this_hit->get_x(1),2)+pow(this_hit->get_y(1),2));
151 
152  lcg += this_hit->get_edep()*(rin + ((rout-rin)/2.0));
153  hcalout_lcg += this_hit->get_edep()*(rin + ((rout-rin)/2.0));
154 
155  }
156 
157  PHG4HitContainer::ConstRange hcalin_abs_hit_range = _hcalin_abs_hit_container->getHits();
158  for (PHG4HitContainer::ConstIterator hit_iter = hcalin_abs_hit_range.first;
159  hit_iter != hcalin_abs_hit_range.second; hit_iter++){
160 
161  PHG4Hit *this_hit = hit_iter->second;
162  assert(this_hit);
163 
164  ea_hcin += this_hit->get_edep();
165 
166  float rin = sqrt(pow(this_hit->get_x(0),2)+pow(this_hit->get_y(0),2));
167  float rout = sqrt(pow(this_hit->get_x(1),2)+pow(this_hit->get_y(1),2));
168 
169  lcg += this_hit->get_edep()*(rin + ((rout-rin)/2.0));
170  hcalin_lcg += this_hit->get_edep()*(rin + ((rout-rin)/2.0));
171 
172  }
173 
174  PHG4HitContainer::ConstRange cemc_abs_hit_range = _cemc_abs_hit_container->getHits();
175  for (PHG4HitContainer::ConstIterator hit_iter = cemc_abs_hit_range.first;
176  hit_iter != cemc_abs_hit_range.second; hit_iter++){
177 
178  PHG4Hit *this_hit = hit_iter->second;
179  assert(this_hit);
180 
181  ea_cemc += this_hit->get_edep();
182 
183  float rin = sqrt(pow(this_hit->get_x(0),2)+pow(this_hit->get_y(0),2));
184  float rout = sqrt(pow(this_hit->get_x(1),2)+pow(this_hit->get_y(1),2));
185 
186  lcg += this_hit->get_edep()*(rin + ((rout-rin)/2.0));
187  emc_lcg += this_hit->get_edep()*(rin + ((rout-rin)/2.0));
188 
189  }
190 
191  PHG4HitContainer::ConstRange magnet_hit_range = _magnet_hit_container->getHits();
192  for (PHG4HitContainer::ConstIterator hit_iter = magnet_hit_range.first;
193  hit_iter != magnet_hit_range.second; hit_iter++){
194 
195  PHG4Hit *this_hit = hit_iter->second;
196  assert(this_hit);
197 
198  e_magnet += this_hit->get_edep();
199 
200  float rin = sqrt(pow(this_hit->get_x(0),2)+pow(this_hit->get_y(0),2));
201  float rout = sqrt(pow(this_hit->get_x(1),2)+pow(this_hit->get_y(1),2));
202 
203  lcg += this_hit->get_edep()*(rin + ((rout-rin)/2.0));
204 
205  }
206 
207  PHG4HitContainer::ConstRange bh_hit_range = _bh_hit_container->getHits();
208  for (PHG4HitContainer::ConstIterator hit_iter = bh_hit_range.first;
209  hit_iter != bh_hit_range.second; hit_iter++){
210 
211  PHG4Hit *this_hit = hit_iter->second;
212  assert(this_hit);
213 
214  e_bh += this_hit->get_edep();
215 
216  }
217 
218  PHG4HitContainer::ConstRange bh_plus_hit_range = _bh_plus_hit_container->getHits();
219  for (PHG4HitContainer::ConstIterator hit_iter = bh_plus_hit_range.first;
220  hit_iter != bh_plus_hit_range.second; hit_iter++){
221 
222  PHG4Hit *this_hit = hit_iter->second;
223  assert(this_hit);
224 
225  e_bh_plus += this_hit->get_edep();
226 
227  }
228 
229  PHG4HitContainer::ConstRange bh_minus_hit_range = _bh_minus_hit_container->getHits();
230  for (PHG4HitContainer::ConstIterator hit_iter = bh_minus_hit_range.first;
231  hit_iter != bh_minus_hit_range.second; hit_iter++){
232 
233  PHG4Hit *this_hit = hit_iter->second;
234  assert(this_hit);
235 
236  e_bh_minus += this_hit->get_edep();
237 
238  }
239 
240  PHG4HitContainer::ConstRange cemc_electronics_hit_range = _cemc_electronics_hit_container->getHits();
241  for (PHG4HitContainer::ConstIterator hit_iter = cemc_electronics_hit_range.first;
242  hit_iter != cemc_electronics_hit_range.second; hit_iter++){
243 
244  PHG4Hit *this_hit = hit_iter->second;
245  assert(this_hit);
246 
247  e_emcelect += this_hit->get_edep();
248 
249  float rin = sqrt(pow(this_hit->get_x(0),2)+pow(this_hit->get_y(0),2));
250  float rout = sqrt(pow(this_hit->get_x(1),2)+pow(this_hit->get_y(1),2));
251 
252  lcg += this_hit->get_edep()*(rin + ((rout-rin)/2.0));
253 
254  }
255 
256  PHG4HitContainer::ConstRange hcalin_spt_hit_range = _hcalin_spt_hit_container->getHits();
257  for (PHG4HitContainer::ConstIterator hit_iter = hcalin_spt_hit_range.first;
258  hit_iter != hcalin_spt_hit_range.second; hit_iter++){
259 
260  PHG4Hit *this_hit = hit_iter->second;
261  assert(this_hit);
262 
263  e_hcalin_spt += this_hit->get_edep();
264 
265  float rin = sqrt(pow(this_hit->get_x(0),2)+pow(this_hit->get_y(0),2));
266  float rout = sqrt(pow(this_hit->get_x(1),2)+pow(this_hit->get_y(1),2));
267 
268  lcg += this_hit->get_edep()*(rin + ((rout-rin)/2.0));
269 
270  }
271 
272  // SVTX
273 
274  PHG4HitContainer::ConstRange svtx_hit_range = _svtx_hit_container->getHits();
275  for (PHG4HitContainer::ConstIterator hit_iter = svtx_hit_range.first;
276  hit_iter != svtx_hit_range.second; hit_iter++){
277 
278  PHG4Hit *this_hit = hit_iter->second;
279  assert(this_hit);
280 
281  e_svtx += this_hit->get_edep();
282 
283  }
284 
285  PHG4HitContainer::ConstRange svtx_support_hit_range = _svtx_support_hit_container->getHits();
286  for (PHG4HitContainer::ConstIterator hit_iter = svtx_support_hit_range.first;
287  hit_iter != svtx_support_hit_range.second; hit_iter++){
288 
289  PHG4Hit *this_hit = hit_iter->second;
290  assert(this_hit);
291 
292  e_svtx_support += this_hit->get_edep();
293 
294  }
295 
296  //normalize the LCG
297 
298  lcg = lcg/(e_cemc+ea_cemc+e_emcelect+e_hcin+ea_hcin+e_hcalin_spt+e_hcout+ea_hcout+e_magnet);
299  hcalin_lcg = hcalin_lcg/(e_hcin+ea_hcin);
300  hcalout_lcg = hcalout_lcg/(e_hcout+ea_hcout);
301  emc_lcg = emc_lcg/(e_cemc+ea_cemc);
302 
303  // Fill the ntuple
304 
305  float ntdata[27];
306 
307  ntdata[0] = e_hcin;
308  ntdata[1] = e_hcout;
309  ntdata[2] = e_cemc;
310  ntdata[3] = ea_hcin;
311  ntdata[4] = ea_hcout;
312  ntdata[5] = ea_cemc;
313  ntdata[6] = ev_hcin;
314  ntdata[7] = ev_hcout;
315  ntdata[8] = ev_cemc;
316  ntdata[9] = e_magnet;
317  ntdata[10] = e_bh;
318  ntdata[11] = e_emcelect;
319  ntdata[12] = e_hcalin_spt;
320 
321  ntdata[13] = e_cemc/(e_cemc + ea_cemc + e_emcelect);
322  ntdata[14] = e_hcin/(e_hcin + ea_hcin + e_hcalin_spt);
323  ntdata[15] = e_hcout/(e_hcout + ea_hcout);
324 
325  ntdata[16] = ev_cemc/(e_cemc + ea_cemc + e_emcelect);
326  ntdata[17] = ev_hcin/(e_hcin + ea_hcin + e_hcalin_spt);
327  ntdata[18] = ev_hcout/(e_hcout + ea_hcout);
328 
329  ntdata[19] = lcg;
330 
331  ntdata[20] = emc_lcg;
332  ntdata[21] = hcalin_lcg;
333  ntdata[22] = hcalout_lcg;
334 
335  ntdata[23] = e_svtx;
336  ntdata[24] = e_svtx_support;
337 
338  ntdata[25] = e_bh_plus;
339  ntdata[26] = e_bh_minus;
340 
341  calenergy->Fill(ntdata);
342 
343  // crosscheck
344 
345  //std::cout << std::endl;
346  //std::cout << "ev_cemc = " << ev_cemc << " e_cemc + ea_cemc + e_emcelect = " << e_cemc + ea_cemc + e_emcelect << " ev_cemc/sfvemc " << ntdata[8]/ntdata[16] << std::endl;
347  //std::cout << "ev_hcin = " << ev_hcin << " e_hcin + ea_hcin + e_hcalin_spt = " << e_hcin + ea_hcin + e_hcalin_spt << " ev_hcin/sfihcal " << ntdata[7]/ntdata[17] << std::endl;
348  //std::cout << "ev_hcout = " << ev_hcout << " e_hcout + ea_hcout = " << e_hcout + ea_hcout << " ev_hcout/sfohcal " << ntdata[6]/ntdata[18] << std::endl;
349  //std::cout << std::endl;
350 
351 
352  // any other return code might lead to aborting the event or analysis for everyone
353  return 0;
354 }
355 
356 int HCALAnalysis::End(PHCompositeNode *topNode)
357 {
358 
359  outputfile->Write();
360  outputfile->Close();
361 
362  return 0;
363 }
364 
365 void HCALAnalysis::GetNodes(PHCompositeNode *topNode)
366 {
367 
368  // store top node
369  _topNode = topNode;
370 
371  // get DST objects
372  _hcalout_hit_container = findNode::getClass<PHG4HitContainer> (topNode, "G4HIT_HCALOUT");
373  if(!_hcalout_hit_container) { std::cout << PHWHERE << ":: No G4HIT_HCALOUT! No sense continuing" << std::endl; exit(1);}
374 
375  _hcalin_hit_container = findNode::getClass<PHG4HitContainer> (topNode, "G4HIT_HCALIN");
376  if(!_hcalin_hit_container) { std::cout << PHWHERE << ":: No G4HIT_HCALIN! No sense continuing" << std::endl; exit(1);}
377 
378  _cemc_hit_container = findNode::getClass<PHG4HitContainer> (topNode, "G4HIT_CEMC");
379  if(!_cemc_hit_container) { std::cout << PHWHERE << ":: No G4HIT_CEMC! No sense continuing" << std::endl; exit(1);}
380 
381  _hcalout_abs_hit_container = findNode::getClass<PHG4HitContainer> (topNode, "G4HIT_ABSORBER_HCALOUT");
382  if(!_hcalout_abs_hit_container) { std::cout << PHWHERE << ":: No G4HIT_HCALOUT! No sense continuing" << std::endl; exit(1);}
383 
384  _hcalin_abs_hit_container = findNode::getClass<PHG4HitContainer> (topNode, "G4HIT_ABSORBER_HCALIN");
385  if(!_hcalin_abs_hit_container) { std::cout << PHWHERE << ":: No G4HIT_HCALIN! No sense continuing" << std::endl; exit(1);}
386 
387  _cemc_abs_hit_container = findNode::getClass<PHG4HitContainer> (topNode, "G4HIT_ABSORBER_CEMC");
388  if(!_cemc_abs_hit_container) { std::cout << PHWHERE << ":: No G4HIT_CEMC! No sense continuing" << std::endl; exit(1);}
389 
390  _magnet_hit_container = findNode::getClass<PHG4HitContainer> (topNode, "G4HIT_MAGNET");
391  if(!_magnet_hit_container) { std::cout << PHWHERE << ":: No G4HIT_MAGNET! No sense continuing" << std::endl; exit(1);}
392 
393  _bh_hit_container = findNode::getClass<PHG4HitContainer> (topNode, "G4HIT_BH_1");
394  if(!_bh_hit_container) { std::cout << PHWHERE << ":: No G4HIT_BH_1! No sense continuing" << std::endl; exit(1);}
395 
396  _cemc_electronics_hit_container = findNode::getClass<PHG4HitContainer> (topNode, "G4HIT_CEMC_ELECTRONICS");
397  if(!_cemc_electronics_hit_container) { std::cout << PHWHERE << ":: No G4HIT_EMCELECTRONICS_0! No sense continuing" << std::endl; exit(1);}
398 
399  _hcalin_spt_hit_container = findNode::getClass<PHG4HitContainer> (topNode, "G4HIT_HCALIN_SPT");
400  if(!_hcalin_spt_hit_container) { std::cout << PHWHERE << ":: No G4HIT_HCALIN_SPT! No sense continuing" << std::endl; exit(1);}
401 
402  _svtx_hit_container = findNode::getClass<PHG4HitContainer> (topNode, "G4HIT_SVTX");
403  if(!_svtx_hit_container) { std::cout << PHWHERE << ":: No G4HIT_SVTX! No sense continuing" << std::endl; exit(1);}
404 
405  _svtx_support_hit_container = findNode::getClass<PHG4HitContainer> (topNode, "G4HIT_SVTXSUPPORT");
406  if(!_svtx_support_hit_container) { std::cout << PHWHERE << ":: No G4HIT_SVTX_SUPPORT! No sense continuing" << std::endl; exit(1);}
407 
408  _bh_plus_hit_container = findNode::getClass<PHG4HitContainer> (topNode, "G4HIT_BH_FORWARD_PLUS");
409  if(!_bh_plus_hit_container) { std::cout << PHWHERE << ":: No G4HIT_BH_FORWARD_PLUS! No sense continuing" << std::endl; exit(1);}
410 
411  _bh_minus_hit_container = findNode::getClass<PHG4HitContainer> (topNode, "G4HIT_BH_FORWARD_NEG");
412  if(!_bh_minus_hit_container) { std::cout << PHWHERE << ":: No G4HIT_BH_FORWARD_NEG! No sense continuing" << std::endl; exit(1);}
413 
414  return;
415 }
416