Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
EMCalLikelihood.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file EMCalLikelihood.C
1 #include "EMCalLikelihood.h"
2 #include "EMCalTrk.h"
3 
4 #include <phool/getClass.h>
5 #include <fun4all/SubsysReco.h>
6 #include <fun4all/Fun4AllServer.h>
7 #include <fun4all/PHTFileServer.h>
8 #include <phool/PHCompositeNode.h>
9 #include <fun4all/Fun4AllReturnCodes.h>
10 #include <fun4all/Fun4AllHistoManager.h>
11 
12 #include <phool/PHCompositeNode.h>
13 
14 #include <g4main/PHG4TruthInfoContainer.h>
15 #include <g4main/PHG4VtxPoint.h>
16 #include <g4main/PHG4Particle.h>
17 
18 #include <g4hough/SvtxVertexMap.h>
19 
20 #include <calobase/RawTowerContainer.h>
21 #include <calobase/RawTowerGeom.h>
22 #include <calobase/RawTower.h>
23 #include <calobase/RawClusterContainer.h>
24 #include <calobase/RawCluster.h>
25 
26 #include <g4eval/CaloEvalStack.h>
27 #include <g4eval/CaloRawClusterEval.h>
28 #include <g4eval/CaloRawTowerEval.h>
29 #include <g4eval/CaloTruthEval.h>
30 #include <g4eval/SvtxEvalStack.h>
31 
32 #include <TNtuple.h>
33 #include <TFile.h>
34 #include <TH1F.h>
35 #include <TH2F.h>
36 #include <TVector3.h>
37 #include <TLorentzVector.h>
38 
39 #include <exception>
40 #include <stdexcept>
41 #include <iostream>
42 #include <vector>
43 #include <set>
44 #include <algorithm>
45 #include <cassert>
46 #include <cmath>
47 
48 using namespace std;
49 
50 EMCalLikelihood::EMCalLikelihood(const std::string &filename) :
51  SubsysReco("EMCalLikelihood"), _filename(filename), _ievent(0), _trk(NULL) //
52  , center_cemc_iphi(NAN), center_cemc_ieta(NAN), //
53  center_hcalin_iphi(NAN), center_hcalin_ieta(NAN), //
54  width_cemc_iphi(NAN), width_cemc_ieta(NAN), //
55  width_hcalin_iphi(NAN), width_hcalin_ieta(NAN), //
56  h2_Edep_Distribution_e(NULL), h2_Edep_Distribution_pi(NULL), //
57  h1_ep_Distribution_e(NULL), h1_ep_Distribution_pi(NULL), //
58  _do_ganging(false), _ganging_size(1, 1)
59 {
60 
61 }
62 
63 EMCalLikelihood::~EMCalLikelihood()
64 {
65 
66 }
67 
68 int
69 EMCalLikelihood::InitRun(PHCompositeNode *topNode)
70 {
71  _ievent = 0;
72 
73  PHNodeIterator iter(topNode);
74  PHCompositeNode *dstNode = static_cast<PHCompositeNode*>(iter.findFirst(
75  "PHCompositeNode", "DST"));
76  if (!dstNode)
77  {
78  std::cerr << PHWHERE << "DST Node missing, doing nothing." << std::endl;
79  throw runtime_error(
80  "Failed to find DST node in EmcRawTowerBuilder::CreateNodes");
81  }
82 
83  {
84  _trk = findNode::getClass<EMCalTrk>(dstNode, "EMCalTrk");
85 
86  if (not _trk)
87  {
88 
89  cout << __PRETTY_FUNCTION__ << "::" << Name() << " - Error - "
90  << "Can not find the EMCalTrk node" << endl;
91 
92  return Fun4AllReturnCodes::ABORTRUN;
93 
94  }
95 
96  assert(_trk);
97  }
98 
99  return Fun4AllReturnCodes::EVENT_OK;
100 }
101 
102 int
103 EMCalLikelihood::End(PHCompositeNode *topNode)
104 {
105 
106  cout << "EMCalLikelihood::End - write to " << _filename << endl;
107  PHTFileServer::get().cd(_filename);
108 
109  Fun4AllHistoManager *hm = get_HistoManager();
110  assert(hm);
111  for (unsigned int i = 0; i < hm->nHistos(); i++)
112  hm->getHisto(i)->Write();
113 
114  return Fun4AllReturnCodes::EVENT_OK;
115 }
116 
117 int
118 EMCalLikelihood::Init(PHCompositeNode *topNode)
119 {
120 
121  cout << "EMCalLikelihood::Init - Making PHTFileServer " << _filename << endl;
122  PHTFileServer::get().open(_filename, "RECREATE");
123 
124  Fun4AllHistoManager *hm = get_HistoManager();
125  assert(hm);
126 
127  TH2F * h_merge_direction = new TH2F("h_merge_direction",
128  "h_merge_direction;merge eta shift;merge phi shift", 10, -.5, 9.5, 10,
129  -.5, 9.5);
130  hm->registerHisto(h_merge_direction);
131 
132  if (h2_Edep_Distribution_e)
133  {
134  h2_Edep_Distribution_e =
135  dynamic_cast<TH2 *>(h2_Edep_Distribution_e->Clone(
136  "h2_Edep_Distribution_e"));
137  hm->registerHisto(h2_Edep_Distribution_e);
138 
139  h1_ep_Distribution_e = h2_Edep_Distribution_e->ProjectionX(
140  "h1_ep_Distribution_e");
141  assert(h1_ep_Distribution_e);
142  hm->registerHisto(h1_ep_Distribution_e);
143 
144  // regulate minimal probability value
145  double min_prob = 1;
146  for (int x = 1; x <= h2_Edep_Distribution_e->GetNbinsX(); ++x)
147  for (int y = 1; y <= h2_Edep_Distribution_e->GetNbinsX(); ++y)
148  {
149  const double binc = h2_Edep_Distribution_e->GetBinContent(x, y);
150  assert(binc >= 0);
151 
152  if (binc > 0 and binc < min_prob)
153  min_prob = binc;
154  }
155  assert(min_prob > 0);
156 
157  for (int x = 0; x <= h2_Edep_Distribution_e->GetNbinsX() + 1; ++x)
158  for (int y = 0; y <= h2_Edep_Distribution_e->GetNbinsX() + 1; ++y)
159  {
160  const double binc = h2_Edep_Distribution_e->GetBinContent(x, y);
161 
162  if (binc == 0)
163  h2_Edep_Distribution_e->SetBinContent(x, y, min_prob);
164  }
165 
166  for (int x = 0; x <= h1_ep_Distribution_e->GetNbinsX() + 1; ++x)
167  {
168  const double binc = h1_ep_Distribution_e->GetBinContent(x);
169 
170  if (binc == 0)
171  h1_ep_Distribution_e->SetBinContent(x, min_prob);
172  }
173  }
174  if (h2_Edep_Distribution_pi)
175  {
176  h2_Edep_Distribution_pi =
177  dynamic_cast<TH2 *>(h2_Edep_Distribution_pi->Clone(
178  "h2_Edep_Distribution_pi"));
179  hm->registerHisto(h2_Edep_Distribution_pi);
180 
181  h1_ep_Distribution_pi = h2_Edep_Distribution_pi->ProjectionX(
182  "h1_ep_Distribution_pi");
183  assert(h1_ep_Distribution_pi);
184  hm->registerHisto(h1_ep_Distribution_pi);
185 
186  // regulate minimal probability value
187  double min_prob = 1;
188  for (int x = 1; x <= h2_Edep_Distribution_pi->GetNbinsX(); ++x)
189  for (int y = 1; y <= h2_Edep_Distribution_pi->GetNbinsX(); ++y)
190  {
191  const double binc = h2_Edep_Distribution_pi->GetBinContent(x, y);
192  assert(binc >= 0);
193 
194  if (binc > 0 and binc < min_prob)
195  min_prob = binc;
196  }
197  assert(min_prob > 0);
198 
199  for (int x = 0; x <= h2_Edep_Distribution_pi->GetNbinsX() + 1; ++x)
200  for (int y = 0; y <= h2_Edep_Distribution_pi->GetNbinsX() + 1; ++y)
201  {
202  const double binc = h2_Edep_Distribution_pi->GetBinContent(x, y);
203 
204  if (binc == 0)
205  h2_Edep_Distribution_pi->SetBinContent(x, y, min_prob);
206  }
207 
208  for (int x = 0; x <= h1_ep_Distribution_pi->GetNbinsX() + 1; ++x)
209  {
210  const double binc = h1_ep_Distribution_pi->GetBinContent(x);
211 
212  if (binc == 0)
213  h1_ep_Distribution_pi->SetBinContent(x, min_prob);
214  }
215  }
216  return Fun4AllReturnCodes::EVENT_OK;
217 }
218 
219 int
220 EMCalLikelihood::process_event(PHCompositeNode *topNode)
221 {
222 // const double significand = _ievent / TMath::Power(10, (int) (log10(_ievent)));
223 
224 // if (fmod(significand, 1.0) == 0 && significand <= 10)
225 // cout << "EMCalLikelihood::process_event - " << _ievent << endl;
226 
227  _ievent++;
228 
229  if (_do_ganging)
230  ApplyEMCalGanging(_trk);
231 
232  UpdateEnergyDeposition(_trk);
233 
234  if (h2_Edep_Distribution_e and h2_Edep_Distribution_pi)
235  UpdateEnergyDepositionLikelihood(_trk);
236 
237  return Fun4AllReturnCodes::EVENT_OK;
238 }
239 
240 Fun4AllHistoManager *
241 EMCalLikelihood::get_HistoManager()
242 {
243 
244  Fun4AllServer *se = Fun4AllServer::instance();
245  Fun4AllHistoManager *hm = se->getHistoManager("EMCalLikelihood_HISTOS");
246 
247  if (not hm)
248  {
249  cout
250  << "EMCalLikelihood::get_HistoManager - Making Fun4AllHistoManager EMCalLikelihood_HISTOS"
251  << endl;
252  hm = new Fun4AllHistoManager("EMCalLikelihood_HISTOS");
253  se->registerHistoManager(hm);
254  }
255 
256  assert(hm);
257 
258  return hm;
259 }
260 
261 void
262 EMCalLikelihood::ApplyEMCalGanging(EMCalTrk * trk)
263 {
264  assert(trk);
265  assert(_ganging_size.first > 1 or _ganging_size.second > 1);
266  assert(_ganging_size.first > 0 and _ganging_size.second > 0);
267 
268  static bool once = true;
269  if (once)
270  {
271  once = false;
272 
273  cout << "EMCalLikelihood::ApplyEMCalGanging - apply EMCal ganging "
274  << _ganging_size.first << "x" << _ganging_size.second << endl;
275  }
276 
277  // estimate an eta-phi bin
278  TVector3 vertex(trk->gvx, trk->gvy, trk->gvz);
279  TVector3 momentum(trk->gpx, trk->gpy, trk->gpz);
280  const double radius = 100; // center of CEMC.
281  const double eta_bin_cm = 2.5;
282  const double phi_bin_rad = 0.025;
283 
284  assert(momentum.Pt() > 0);
285 
286  TVector3 proj = vertex + radius / momentum.Pt() * momentum;
287  const int eta_bin = floor(proj.Z() / eta_bin_cm);
288  const int phi_bin = floor(proj.Phi() / phi_bin_rad);
289 
290  const int merge_plus_eta = (eta_bin % _ganging_size.first);
291  const bool merge_plus_phi = (phi_bin % _ganging_size.second);
292 
293  Fun4AllHistoManager *hm = get_HistoManager();
294  assert(hm);
295 
296  TH2F * h_merge_direction = (TH2F *) hm->getHisto("h_merge_direction");
297  assert(h_merge_direction);
298 
299  h_merge_direction->Fill(merge_plus_eta, merge_plus_phi);
300 
301  for (int ieta = 0; ieta < trk->Max_N_Tower; ++ieta)
302  {
303  for (int iphi = 0; iphi < trk->Max_N_Tower; ++iphi)
304  {
305  bool out_of_edge = false;
306 
307  double cemc_ieta = 0;
308  double cemc_iphi = 0;
309  double cemc_energy = 0;
310 
311  for (unsigned int dieta = 0; dieta < _ganging_size.first; ++dieta)
312  {
313  const unsigned int fetch_eta = _ganging_size.first * ieta + dieta + merge_plus_eta;
314 
315  if ( fetch_eta>= trk->Max_N_Tower or out_of_edge)
316  {
317  out_of_edge = true;
318  break;
319  }
320 
321  for (unsigned int diphi = 0; diphi < _ganging_size.second; ++diphi)
322  {
323  const unsigned int fetch_phi = _ganging_size.first * iphi + diphi + merge_plus_phi;
324 
325  if (fetch_phi>= trk->Max_N_Tower or out_of_edge)
326  {
327  out_of_edge = true;
328  break;
329  }
330 
331  cemc_ieta += trk->cemc_ieta[fetch_eta][fetch_phi];
332  cemc_iphi += trk->cemc_iphi[fetch_eta][fetch_phi];
333  cemc_energy += trk->cemc_energy[fetch_eta][fetch_phi];
334  }
335  }
336 
337  if (out_of_edge)
338  {
339  trk->cemc_ieta[ieta][iphi] = -9999;
340  trk->cemc_iphi[ieta][iphi] = -9999;
341  trk->cemc_energy[ieta][iphi] = -0;
342  }
343  else
344  {
345  trk->cemc_ieta[ieta][iphi] = cemc_ieta / _ganging_size.first
346  / _ganging_size.second;
347  trk->cemc_iphi[ieta][iphi] = cemc_iphi / _ganging_size.first
348  / _ganging_size.second;
349  trk->cemc_energy[ieta][iphi] = cemc_energy;
350  }
351  }
352  }
353 }
354 
355 void
356 EMCalLikelihood::UpdateEnergyDeposition(EMCalTrk * trk)
357 {
358  assert(trk);
359 
360  trk->cemc_sum_energy = 0;
361  trk->cemc_sum_n_tower = 0;
362  trk->hcalin_sum_energy = 0;
363  trk->hcalin_sum_n_tower = 0;
364 
365  for (int ieta = 0; ieta < trk->Max_N_Tower; ++ieta)
366  {
367  for (int iphi = 0; iphi < trk->Max_N_Tower; ++iphi)
368  {
369 
370  //cemc
371  const double cemc_radius2 = pow(
372  (trk->cemc_ieta[ieta][iphi] - center_cemc_ieta) / width_cemc_ieta,
373  2.)
374  + pow(
375  (trk->cemc_iphi[ieta][iphi] - center_cemc_iphi)
376  / width_cemc_iphi, 2.);
377  trk->cemc_radius[ieta][iphi] = sqrt(cemc_radius2);
378 
379  if (cemc_radius2 <= 1.0)
380  {
381  trk->cemc_sum_energy += trk->cemc_energy[ieta][iphi];
382  trk->cemc_sum_n_tower++;
383  }
384 
385  // hcalin
386  const double hcalin_radius2 = pow(
387  (trk->hcalin_ieta[ieta][iphi] - center_hcalin_ieta)
388  / width_hcalin_ieta, 2.)
389  + pow(
390  (trk->hcalin_iphi[ieta][iphi] - center_hcalin_iphi)
391  / width_hcalin_iphi, 2.);
392  trk->hcalin_radius[ieta][iphi] = sqrt(hcalin_radius2);
393 
394  if (hcalin_radius2 <= 1.0)
395  {
396  trk->hcalin_sum_energy += trk->hcalin_energy[ieta][iphi];
397  trk->hcalin_sum_n_tower++;
398  }
399  }
400  }
401 }
402 
403 void
404 EMCalLikelihood::UpdateEnergyDepositionLikelihood(EMCalTrk * trk)
405 {
406  assert(trk);
407  assert(h2_Edep_Distribution_e);
408  assert(h2_Edep_Distribution_pi);
409  assert(h1_ep_Distribution_e);
410  assert(h1_ep_Distribution_pi);
411 
412  assert(
413  h2_Edep_Distribution_e->GetNbinsX() == h2_Edep_Distribution_pi->GetNbinsX());
414  assert(
415  h2_Edep_Distribution_e->GetNbinsY() == h2_Edep_Distribution_pi->GetNbinsY());
416  assert(
417  h2_Edep_Distribution_e->GetNbinsX() == h1_ep_Distribution_e->GetNbinsX());
418  assert(
419  h1_ep_Distribution_pi->GetNbinsX() == h2_Edep_Distribution_pi->GetNbinsX());
420 
421  const int binx = h2_Edep_Distribution_e->GetXaxis()->FindBin(trk->get_ep());
422  const int biny = h2_Edep_Distribution_e->GetYaxis()->FindBin(
423  trk->hcalin_sum_energy);
424 
425  {
426  double prob_e = h2_Edep_Distribution_e->GetBinContent(binx, biny);
427  if (!prob_e > 0)
428  {
429  cout << __PRETTY_FUNCTION__ << Name()
430  << " - Error - invalid likelihood value prob_e = " << prob_e
431  << " @ bin " << binx << ", " << biny << endl;
432  }
433  assert(prob_e > 0);
434 
435  double prob_pi = h2_Edep_Distribution_pi->GetBinContent(binx, biny);
436  if (!prob_pi > 0)
437  {
438  cout << __PRETTY_FUNCTION__ << Name()
439  << " - Error - invalid likelihood value prob_pi = " << prob_pi
440  << " @ bin " << binx << ", " << biny << endl;
441  }
442  assert(prob_pi > 0);
443 
444  trk->ll_edep_e = log(prob_e);
445  trk->ll_edep_h = log(prob_pi);
446  }
447 
448  {
449  double prob_e = h1_ep_Distribution_e->GetBinContent(binx);
450  assert(prob_e > 0);
451 
452  double prob_pi = h1_ep_Distribution_pi->GetBinContent(binx);
453  assert(prob_pi > 0);
454 
455  trk->ll_ep_e = log(prob_e);
456  trk->ll_ep_h = log(prob_pi);
457  }
458 }