Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
ShowerSize.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file ShowerSize.cc
1 #include "ShowerSize.h"
2 
3 #include <g4main/PHG4HitContainer.h>
4 #include <g4main/PHG4Hit.h>
5 
6 #include <g4main/PHG4TruthInfoContainer.h>
7 #include <g4main/PHG4Particle.h>
8 #include <g4main/PHG4VtxPoint.h>
9 
10 #include <fun4all/Fun4AllHistoManager.h>
11 
12 #include <phool/getClass.h>
13 
14 #include <TFile.h>
15 #include <TH1.h>
16 #include <TH2.h>
17 #include <TNtuple.h>
18 
19 #include<sstream>
20 
21 using namespace std;
22 
23 ShowerSize::ShowerSize(const std::string &name, const std::string &filename):
24  SubsysReco( name ),
25  nblocks(0),
26  hm(nullptr),
27  _filename(filename),
28  ntups(nullptr),
29  ntupe(nullptr),
30  ntup(nullptr),
31  outfile(nullptr)
32 {}
33 
34 ShowerSize::~ShowerSize()
35 {
36  // delete ntup;
37  delete hm;
38 }
39 
40 
41 int
42 ShowerSize::Init( PHCompositeNode* )
43 {
44  ostringstream hname, htit;
45  hm = new Fun4AllHistoManager(Name());
46  outfile = new TFile(_filename.c_str(), "RECREATE");
47  ntups = new TNtuple("sz"," Shower size","rad:em:hi:ho:mag:bh");
48  ntupe = new TNtuple("truth", "The Absolute Truth", "phi:theta:eta:e:p");
49  ntup = new TNtuple("de", "Change in Angles", "ID:dphi:dtheta:dtotal:edep");
50  return 0;
51 }
52 
53 int
54 ShowerSize::process_event( PHCompositeNode* topNode )
55 {
56  map<int, PHG4Particle*>::const_iterator particle_iter;
57  PHG4TruthInfoContainer *_truth_container = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
58 
59  PHG4TruthInfoContainer::ConstRange primary_range =
60  _truth_container->GetPrimaryParticleRange();
61  float ntvars[5] = {0};
62  for (PHG4TruthInfoContainer::ConstIterator particle_iter = primary_range.first;
63  particle_iter != primary_range.second; ++particle_iter)
64  {
65  PHG4Particle *particle = particle_iter->second;
66  ntvars[0] = atan2(particle->get_py(), particle->get_px());
67  ntvars[1] = atan(sqrt(particle->get_py()*particle->get_py() +
68  particle->get_px()*particle->get_px()) /
69  particle->get_pz());
70  if (ntvars[1] < 0)
71  {
72  ntvars[1]+=M_PI;
73  }
74  ntvars[2] = 0.5*log((particle->get_e()+particle->get_pz())/
75  (particle->get_e()-particle->get_pz()));
76  ntvars[3] = particle->get_e();
77  ntvars[4] = particle->get_pz();
78  }
79  ntupe->Fill(ntvars);
80  PHG4HitContainer *hits = nullptr;
81  double emc[10] = {0};
82  double ih[10] = {0};
83  double oh[10] = {0};
84  double mag[10] = {0};
85  double bh[10] = {0};
86  double eall[5] = {0};
87  float detid[5] = {0};
88  string nodename[2] = {"G4HIT_CEMC","G4HIT_ABSORBER_CEMC"};
89  for (int j=0; j<2;j++)
90  {
91  hits = findNode::getClass<PHG4HitContainer>(topNode,nodename[j]);
92  if (hits)
93  {
94  PHG4HitContainer::ConstRange hit_range = hits->getHits();
95  for ( PHG4HitContainer::ConstIterator hit_iter = hit_range.first ; hit_iter != hit_range.second; hit_iter++ )
96  {
97  // double radius = sqrt(hit_iter->second->get_avg_x() * hit_iter->second->get_avg_x() +
98  // hit_iter->second->get_avg_y() * hit_iter->second->get_avg_y() +
99  // hit_iter->second->get_avg_z() * hit_iter->second->get_avg_z());
100  double phi = atan2(hit_iter->second->get_avg_y(),hit_iter->second->get_avg_x());
101  double theta = atan(sqrt(hit_iter->second->get_avg_x() * hit_iter->second->get_avg_x() +
102  hit_iter->second->get_avg_y() * hit_iter->second->get_avg_y()) /
103  hit_iter->second->get_avg_z());
104 
105 // handle rollover from pi to -pi
106  double diffphi = phi-ntvars[0];
107  if (diffphi > M_PI)
108  {
109  diffphi -= 2*M_PI;
110  }
111  else if (diffphi < - M_PI)
112  {
113  diffphi += 2*M_PI;
114  }
115  // double difftheta = theta-ntvars[1];
116  if (theta < 0)
117  {
118  theta += M_PI;
119  }
120  double difftheta = theta-ntvars[1];
121 // theta goes from 0-PI --> no rollover problem
122  double deltasqrt = sqrt(diffphi*diffphi+difftheta*difftheta);
123  double edep = hit_iter->second->get_edep();
124  eall[0] += edep;
125  detid[0] = 0;
126  detid[1] = diffphi;
127  detid[2] = theta-ntvars[1];
128  detid[3] = deltasqrt;
129  detid[4] = edep;
130  ntup->Fill(detid);
131  for (int i=0; i<10; i++)
132  {
133  if (deltasqrt < (i+1)*0.025)
134  {
135  emc[i]+=edep;
136  }
137  }
138  }
139  }
140  }
141  nodename[0] = {"G4HIT_HCALIN"};
142  nodename[1] = {"G4HIT_ABSORBER_HCALIN"};
143  for (int j=0; j<2;j++)
144  {
145  hits = findNode::getClass<PHG4HitContainer>(topNode,nodename[j]);
146  if (hits)
147  {
148  PHG4HitContainer::ConstRange hit_range = hits->getHits();
149  for ( PHG4HitContainer::ConstIterator hit_iter = hit_range.first ; hit_iter != hit_range.second; hit_iter++ )
150  {
151  // double radius = sqrt(hit_iter->second->get_avg_x() * hit_iter->second->get_avg_x() +
152  // hit_iter->second->get_avg_y() * hit_iter->second->get_avg_y() +
153  // hit_iter->second->get_avg_z() * hit_iter->second->get_avg_z());
154  double phi = atan2(hit_iter->second->get_avg_y(),hit_iter->second->get_avg_x());
155  double theta = atan(sqrt(hit_iter->second->get_avg_x() * hit_iter->second->get_avg_x() +
156  hit_iter->second->get_avg_y() * hit_iter->second->get_avg_y()) /
157  hit_iter->second->get_avg_z());
158  // handle rollover from pi to -pi
159  double diffphi = phi-ntvars[0];
160  if (diffphi > M_PI)
161  {
162  diffphi -= 2*M_PI;
163  }
164  else if (diffphi < - M_PI)
165  {
166  diffphi += 2*M_PI;
167  }
168  // double difftheta = theta-ntvars[1];
169  if (theta < 0)
170  {
171  theta += M_PI;
172  }
173  double difftheta = theta-ntvars[1];
174 // theta goes from 0-PI --> no rollover problem
175  double deltasqrt = sqrt(diffphi*diffphi+difftheta*difftheta);
176  double edep = hit_iter->second->get_edep();
177  eall[0] += edep;
178  detid[0] = 1;
179  detid[1] = diffphi;
180  detid[2] = theta-ntvars[1];
181  detid[3] = deltasqrt;
182  detid[4] = edep;
183  ntup->Fill(detid);
184  for (int i=0; i<10; i++)
185  {
186  if (deltasqrt < (i+1)*0.025)
187  {
188  ih[i]+=edep;
189  }
190  }
191  }
192  }
193  }
194  nodename[0] = {"G4HIT_HCALOUT"};
195  nodename[1] = {"G4HIT_ABSORBER_HCALOUT"};
196  for (int j=0; j<2;j++)
197  {
198  hits = findNode::getClass<PHG4HitContainer>(topNode,nodename[j]);
199  if (hits)
200  {
201  PHG4HitContainer::ConstRange hit_range = hits->getHits();
202  for ( PHG4HitContainer::ConstIterator hit_iter = hit_range.first ; hit_iter != hit_range.second; hit_iter++ )
203  {
204  // double radius = sqrt(hit_iter->second->get_avg_x() * hit_iter->second->get_avg_x() +
205  // hit_iter->second->get_avg_y() * hit_iter->second->get_avg_y() +
206  // hit_iter->second->get_avg_z() * hit_iter->second->get_avg_z());
207  double phi = atan2(hit_iter->second->get_avg_y(),hit_iter->second->get_avg_x());
208  double theta = atan(sqrt(hit_iter->second->get_avg_x() * hit_iter->second->get_avg_x() +
209  hit_iter->second->get_avg_y() * hit_iter->second->get_avg_y()) /
210  hit_iter->second->get_avg_z());
211 // handle rollover from pi to -pi
212  double diffphi = phi-ntvars[0];
213  if (diffphi > M_PI)
214  {
215  diffphi -= 2*M_PI;
216  }
217  else if (diffphi < - M_PI)
218  {
219  diffphi += 2*M_PI;
220  }
221  // double difftheta = theta-ntvars[1];
222  if (theta < 0)
223  {
224  theta += M_PI;
225  }
226  double difftheta = theta-ntvars[1];
227 // theta goes from 0-PI --> no rollover problem
228  double deltasqrt = sqrt(diffphi*diffphi+difftheta*difftheta);
229  double edep = hit_iter->second->get_edep();
230  eall[0] += edep;
231  detid[0] = 2;
232  detid[1] = diffphi;
233  detid[2] = theta-ntvars[1];
234  detid[3] = deltasqrt;
235  detid[4] = edep;
236  ntup->Fill(detid);
237  for (int i=0; i<10; i++)
238  {
239  if (deltasqrt < (i+1)*0.025)
240  {
241  oh[i]+=edep;
242  }
243  }
244  }
245  }
246  }
247  hits = findNode::getClass<PHG4HitContainer>(topNode,"G4HIT_MAGNET");
248  if (hits)
249  {
250  PHG4HitContainer::ConstRange hit_range = hits->getHits();
251  for ( PHG4HitContainer::ConstIterator hit_iter = hit_range.first ; hit_iter != hit_range.second; hit_iter++ )
252  {
253  // double radius = sqrt(hit_iter->second->get_avg_x() * hit_iter->second->get_avg_x() +
254  // hit_iter->second->get_avg_y() * hit_iter->second->get_avg_y() +
255  // hit_iter->second->get_avg_z() * hit_iter->second->get_avg_z());
256  double phi = atan2(hit_iter->second->get_avg_y(),hit_iter->second->get_avg_x());
257  double theta = atan(sqrt(hit_iter->second->get_avg_x() * hit_iter->second->get_avg_x() +
258  hit_iter->second->get_avg_y() * hit_iter->second->get_avg_y()) /
259  hit_iter->second->get_avg_z());
260 // handle rollover from pi to -pi
261  double diffphi = phi-ntvars[0];
262  if (diffphi > M_PI)
263  {
264  diffphi -= 2*M_PI;
265  }
266  else if (diffphi < - M_PI)
267  {
268  diffphi += 2*M_PI;
269  }
270  //double difftheta = theta-ntvars[1];
271  if (theta < 0)
272  {
273  theta += M_PI;
274  }
275  double difftheta = theta-ntvars[1];
276 // theta goes from 0-PI --> no rollover problem
277  double deltasqrt = sqrt(diffphi*diffphi+difftheta*difftheta);
278  double edep = hit_iter->second->get_edep();
279  eall[0] += edep;
280  detid[0] = 3;
281  detid[1] = diffphi;
282  detid[2] = theta-ntvars[1];
283  detid[3] = deltasqrt;
284  detid[4] = edep;
285  ntup->Fill(detid);
286  for (int i=0; i<10; i++)
287  {
288  if (deltasqrt < (i+1)*0.025)
289  {
290  mag[i]+=edep;
291  }
292  }
293  }
294  }
295  hits = findNode::getClass<PHG4HitContainer>(topNode,"G4HIT_BH_1");
296  if (hits)
297  {
298  PHG4HitContainer::ConstRange hit_range = hits->getHits();
299  for ( PHG4HitContainer::ConstIterator hit_iter = hit_range.first ; hit_iter != hit_range.second; hit_iter++ )
300  {
301  // double radius = sqrt(hit_iter->second->get_avg_x() * hit_iter->second->get_avg_x() +
302  // hit_iter->second->get_avg_y() * hit_iter->second->get_avg_y() +
303  // hit_iter->second->get_avg_z() * hit_iter->second->get_avg_z());
304  double phi = atan2(hit_iter->second->get_avg_y(),hit_iter->second->get_avg_x());
305  double theta = atan(sqrt(hit_iter->second->get_avg_x() * hit_iter->second->get_avg_x() +
306  hit_iter->second->get_avg_y() * hit_iter->second->get_avg_y()) /
307  hit_iter->second->get_avg_z());
308 // handle rollover from pi to -pi
309  double diffphi = phi-ntvars[0];
310  if (diffphi > M_PI)
311  {
312  diffphi -= 2*M_PI;
313  }
314  else if (diffphi < - M_PI)
315  {
316  diffphi += 2*M_PI;
317  }
318  // double difftheta = theta-ntvars[1];
319  if (theta < 0)
320  {
321  theta += M_PI;
322  }
323  double difftheta = theta-ntvars[1];
324 // theta goes from 0-PI --> no rollover problem
325  double deltasqrt = sqrt(diffphi*diffphi+difftheta*difftheta);
326  double edep = hit_iter->second->get_edep();
327  eall[0] += edep;
328  detid[0] = 4;
329  detid[1] = diffphi;
330  detid[2] = theta-ntvars[1];
331  detid[3] = deltasqrt;
332  detid[4] = edep;
333  ntup->Fill(detid);
334  for (int i=0; i<10; i++)
335  {
336  if (deltasqrt < (i+1)*0.025)
337  {
338  bh[i]+=edep;
339  }
340  }
341  }
342  }
343  for (int i=0; i<10;i++)
344  {
345  float nte[6] = {0};
346  nte[0] = i;
347  nte[1] = emc[i];
348  nte[2] = ih[i];
349  nte[3] = oh[i];
350  nte[4] = mag[i];
351  nte[5] = bh[i];
352  ntups->Fill(nte);
353  }
354  float nte[6] = {0};
355  nte[0] = 100;
356  nte[1] = eall[0];
357  nte[2] = eall[1];
358  nte[3] = eall[2];
359  nte[4] = eall[3];
360  nte[5] = eall[4];
361  ntups->Fill(nte);
362 
363  return 0;
364 }
365 
366 int
367 ShowerSize::End(PHCompositeNode * topNode)
368 {
369  outfile->cd();
370  // ntup->Write();
371  ntups->Write();
372  ntupe->Write();
373  outfile->Write();
374  outfile->Close();
375  delete outfile;
376  hm->dumpHistos(_filename, "UPDATE");
377  return 0;
378 }
379 
380 void
381 ShowerSize::AddNode(const std::string &name, const int detid)
382 {
383  _node_postfix.insert(name);
384  _detid[name] = detid;
385  return;
386 }