Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
EMCalAna.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file EMCalAna.C
1 #include "EMCalAna.h"
2 #include "UpsilonPair.h"
3 
4 #include <fun4all/Fun4AllHistoManager.h>
5 #include <fun4all/Fun4AllReturnCodes.h>
6 #include <fun4all/Fun4AllServer.h>
7 #include <fun4all/PHTFileServer.h>
8 #include <fun4all/SubsysReco.h>
9 #include <phool/PHCompositeNode.h>
10 #include <phool/getClass.h>
11 
12 #include <phool/PHCompositeNode.h>
13 
14 #include <g4main/PHG4Hit.h>
15 #include <g4main/PHG4HitContainer.h>
16 #include <g4main/PHG4Particle.h>
17 #include <g4main/PHG4TruthInfoContainer.h>
18 #include <g4main/PHG4VtxPoint.h>
19 
20 #include <g4hough/SvtxCluster.h>
21 #include <g4hough/SvtxClusterMap.h>
22 #include <g4hough/SvtxTrack.h>
23 #include <g4hough/SvtxVertexMap.h>
24 
25 #include <g4detectors/PHG4CylinderCellGeomContainer.h>
26 #include <g4detectors/PHG4CylinderCellGeom_Spacalv1.h>
27 #include <g4detectors/PHG4CylinderGeomContainer.h>
28 #include <g4detectors/PHG4CylinderGeom_Spacalv3.h>
29 
30 #include <calobase/RawCluster.h>
31 #include <calobase/RawClusterContainer.h>
32 #include <calobase/RawTower.h>
33 #include <calobase/RawTowerContainer.h>
34 #include <calobase/RawTowerGeomContainer.h>
35 
36 #include <g4eval/CaloEvalStack.h>
37 #include <g4eval/CaloRawClusterEval.h>
38 #include <g4eval/CaloRawTowerEval.h>
39 #include <g4eval/CaloTruthEval.h>
40 #include <g4eval/SvtxEvalStack.h>
41 
42 #include <TFile.h>
43 #include <TH1F.h>
44 #include <TH2F.h>
45 #include <TLorentzVector.h>
46 #include <TNtuple.h>
47 #include <TVector3.h>
48 
49 #include <algorithm>
50 #include <cassert>
51 #include <cmath>
52 #include <exception>
53 #include <iostream>
54 #include <set>
55 #include <stdexcept>
56 #include <vector>
57 
58 using namespace std;
59 
60 EMCalAna::EMCalAna(const std::string &filename, EMCalAna::enu_flags flags)
61  : SubsysReco("EMCalAna")
62  , _eval_stack(NULL)
63  , _T_EMCalTrk(NULL)
64  , _trk(NULL)
65  , //
66  _magfield(1.5)
67  , //
68  _filename(filename)
69  , _flags(flags)
70  , _ievent(0)
71 {
72  verbosity = 1;
73 
74  _hcalout_hit_container = NULL;
75  _hcalin_hit_container = NULL;
76  _cemc_hit_container = NULL;
77  _hcalout_abs_hit_container = NULL;
78  _hcalin_abs_hit_container = NULL;
79  _cemc_abs_hit_container = NULL;
80  _magnet_hit_container = NULL;
81  _bh_hit_container = NULL;
82  _bh_plus_hit_container = NULL;
83  _bh_minus_hit_container = NULL;
84  _cemc_electronics_hit_container = NULL;
85  _hcalin_spt_hit_container = NULL;
86  _svtx_hit_container = NULL;
87  _svtx_support_hit_container = NULL;
88 }
89 
90 EMCalAna::~EMCalAna()
91 {
92  if (_eval_stack)
93  {
94  delete _eval_stack;
95  }
96 }
97 
98 int EMCalAna::InitRun(PHCompositeNode *topNode)
99 {
100  _ievent = 0;
101 
102  PHNodeIterator iter(topNode);
103  PHCompositeNode *dstNode = static_cast<PHCompositeNode *>(iter.findFirst(
104  "PHCompositeNode", "DST"));
105  if (!dstNode)
106  {
107  std::cerr << PHWHERE << "DST Node missing, doing nothing." << std::endl;
108  throw runtime_error(
109  "Failed to find DST node in EmcRawTowerBuilder::CreateNodes");
110  }
111 
112  // get DST objects
113  _hcalout_hit_container = findNode::getClass<PHG4HitContainer>(topNode,
114  "G4HIT_HCALOUT");
115  _hcalin_hit_container = findNode::getClass<PHG4HitContainer>(topNode,
116  "G4HIT_HCALIN");
117 
118  _cemc_hit_container = findNode::getClass<PHG4HitContainer>(topNode,
119  "G4HIT_CEMC");
120 
121  _hcalout_abs_hit_container = findNode::getClass<PHG4HitContainer>(topNode,
122  "G4HIT_ABSORBER_HCALOUT");
123 
124  _hcalin_abs_hit_container = findNode::getClass<PHG4HitContainer>(topNode,
125  "G4HIT_ABSORBER_HCALIN");
126 
127  _cemc_abs_hit_container = findNode::getClass<PHG4HitContainer>(topNode,
128  "G4HIT_ABSORBER_CEMC");
129 
130  if (flag(kProcessUpslisonTrig))
131  {
132  UpsilonPair *upsilon_pair = findNode::getClass<UpsilonPair>(dstNode,
133  "UpsilonPair");
134 
135  if (not upsilon_pair)
136  {
137  upsilon_pair = new UpsilonPair();
138 
139  PHIODataNode<PHObject> *node = new PHIODataNode<PHObject>(
140  upsilon_pair, "UpsilonPair", "PHObject");
141  dstNode->addNode(node);
142  }
143 
144  assert(upsilon_pair);
145  }
146  if (flag(kProcessTrk))
147  {
148  _trk = findNode::getClass<EMCalTrk>(dstNode, "EMCalTrk");
149 
150  if (not _trk)
151  {
152  _trk = new EMCalTrk();
153 
154  PHIODataNode<PHObject> *node = new PHIODataNode<PHObject>(_trk,
155  "EMCalTrk", "PHObject");
156  dstNode->addNode(node);
157  }
158 
159  assert(_trk);
160  }
161 
162  return Fun4AllReturnCodes::EVENT_OK;
163 }
164 
165 int EMCalAna::End(PHCompositeNode *topNode)
166 {
167  cout << "EMCalAna::End - write to " << _filename << endl;
168  PHTFileServer::get().cd(_filename);
169 
170  Fun4AllHistoManager *hm = get_HistoManager();
171  assert(hm);
172  for (unsigned int i = 0; i < hm->nHistos(); i++)
173  hm->getHisto(i)->Write();
174 
175  if (_T_EMCalTrk)
176  _T_EMCalTrk->Write();
177 
178  // help index files with TChain
179  TTree *T_Index = new TTree("T_Index", "T_Index");
180  assert(T_Index);
181  T_Index->Write();
182 
183  return Fun4AllReturnCodes::EVENT_OK;
184 }
185 
186 int EMCalAna::Init(PHCompositeNode *topNode)
187 {
188  _ievent = 0;
189 
190  cout << "EMCalAna::get_HistoManager - Making PHTFileServer " << _filename
191  << endl;
192  PHTFileServer::get().open(_filename, "RECREATE");
193 
194  if (flag(kProcessSF))
195  {
196  cout << "EMCalAna::Init - Process sampling fraction" << endl;
197  Init_SF(topNode);
198  }
199  if (flag(kProcessTower))
200  {
201  cout << "EMCalAna::Init - Process tower occupancies" << endl;
202  Init_Tower(topNode);
203  }
204  if (flag(kProcessTrk))
205  {
206  cout << "EMCalAna::Init - Process trakcs" << endl;
207  Init_Trk(topNode);
208  }
209  if (flag(kProcessUpslisonTrig))
210  {
211  cout << "EMCalAna::Init - Process kProcessUpslisonTrig" << endl;
212  Init_UpslisonTrig(topNode);
213  }
214 
215  return Fun4AllReturnCodes::EVENT_OK;
216 }
217 
218 int EMCalAna::process_event(PHCompositeNode *topNode)
219 {
220  if (verbosity > 2)
221  cout << "EMCalAna::process_event() entered" << endl;
222 
223  if (_eval_stack)
224  _eval_stack->next_event(topNode);
225 
226  if (flag(kProcessUpslisonTrig))
227  {
228  int ret = process_event_UpslisonTrig(topNode);
229  if (ret != Fun4AllReturnCodes::EVENT_OK)
230  return ret;
231  }
232  if (flag(kProcessSF))
233  process_event_SF(topNode);
234  if (flag(kProcessTower))
235  process_event_Tower(topNode);
236  if (flag(kProcessTrk))
237  process_event_Trk(topNode);
238 
239  return Fun4AllReturnCodes::EVENT_OK;
240 }
241 
242 int EMCalAna::Init_UpslisonTrig(PHCompositeNode *topNode)
243 {
244  return Fun4AllReturnCodes::EVENT_OK;
245 }
246 
247 int EMCalAna::process_event_UpslisonTrig(PHCompositeNode *topNode)
248 {
249  if (verbosity > 2)
250  cout << "EMCalAna::process_event_UpslisonTrig() entered" << endl;
251 
252  UpsilonPair *upsilon_pair = findNode::getClass<UpsilonPair>(topNode,
253  "UpsilonPair");
254  assert(upsilon_pair);
255  static const double upsilon_mass = 9460.30e-3;
256 
257  if (!_eval_stack)
258  {
259  _eval_stack = new SvtxEvalStack(topNode);
260  _eval_stack->set_strict(false);
261  }
262  //
263  // SvtxVertexEval* vertexeval = svtxevalstack.get_vertex_eval();
264  // SvtxTrackEval* trackeval = svtxevalstack.get_track_eval();
265  // SvtxClusterEval* clustereval = svtxevalstack.get_cluster_eval();
266  // SvtxHitEval* hiteval = svtxevalstack.get_hit_eval();
267  // SvtxTruthEval* trutheval = svtxevalstack.get_truth_eval();
268  //
269  // SvtxTrackMap* trackmap = findNode::getClass<SvtxTrackMap>(topNode,
270  // "SvtxTrackMap");
271  // assert(trackmap);
272  PHG4TruthInfoContainer *truthinfo =
273  findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
274  assert(truthinfo);
275  PHG4TruthInfoContainer::Map particle_map = truthinfo->GetMap();
276  PHG4TruthInfoContainer::ConstRange primary_range =
277  truthinfo->GetPrimaryParticleRange();
278 
279  set<int> e_pos_candidate;
280  set<int> e_neg_candidate;
281 
282  for (PHG4TruthInfoContainer::ConstIterator iter = primary_range.first;
283  iter != primary_range.second; ++iter)
284  {
285  PHG4Particle *particle = iter->second;
286  assert(particle);
287 
288  if (particle->get_pid() == 11)
289  e_neg_candidate.insert(particle->get_track_id());
290  if (particle->get_pid() == -11)
291  e_pos_candidate.insert(particle->get_track_id());
292  }
293 
294  map<double, pair<int, int> > mass_diff_id_map;
295  for (set<int>::const_iterator i_e_pos = e_pos_candidate.begin();
296  i_e_pos != e_pos_candidate.end(); ++i_e_pos)
297  for (set<int>::const_iterator i_e_neg = e_neg_candidate.begin();
298  i_e_neg != e_neg_candidate.end(); ++i_e_neg)
299  {
300  TLorentzVector vpos(particle_map[*i_e_pos]->get_px(),
301  particle_map[*i_e_pos]->get_py(), particle_map[*i_e_pos]->get_pz(),
302  particle_map[*i_e_pos]->get_e());
303  TLorentzVector vneg(particle_map[*i_e_neg]->get_px(),
304  particle_map[*i_e_neg]->get_py(), particle_map[*i_e_neg]->get_pz(),
305  particle_map[*i_e_neg]->get_e());
306 
307  TLorentzVector invar = vpos + vneg;
308 
309  const double mass = invar.M();
310 
311  const double mass_diff = fabs(mass - upsilon_mass);
312 
313  mass_diff_id_map[mass_diff] = make_pair(*i_e_pos, *i_e_neg);
314  }
315 
316  if (mass_diff_id_map.size() <= 0)
317  return Fun4AllReturnCodes::DISCARDEVENT;
318  else
319  {
320  const pair<int, int> best_upsilon_pair = mass_diff_id_map.begin()->second;
321 
322  //truth mass
323  TLorentzVector gvpos(particle_map[best_upsilon_pair.first]->get_px(),
324  particle_map[best_upsilon_pair.first]->get_py(),
325  particle_map[best_upsilon_pair.first]->get_pz(),
326  particle_map[best_upsilon_pair.first]->get_e());
327  TLorentzVector gvneg(particle_map[best_upsilon_pair.second]->get_px(),
328  particle_map[best_upsilon_pair.second]->get_py(),
329  particle_map[best_upsilon_pair.second]->get_pz(),
330  particle_map[best_upsilon_pair.second]->get_e());
331  TLorentzVector ginvar = gvpos + gvneg;
332  upsilon_pair->gmass = ginvar.M();
333 
334  eval_trk(particle_map[best_upsilon_pair.first], *upsilon_pair->get_trk(0),
335  topNode);
336  eval_trk(particle_map[best_upsilon_pair.second],
337  *upsilon_pair->get_trk(1), topNode);
338 
339  //calculated mass
340  TLorentzVector vpos(upsilon_pair->get_trk(0)->px,
341  upsilon_pair->get_trk(0)->py, upsilon_pair->get_trk(0)->pz, 0);
342  TLorentzVector vneg(upsilon_pair->get_trk(1)->px,
343  upsilon_pair->get_trk(1)->py, upsilon_pair->get_trk(1)->pz, 0);
344  vpos.SetE(vpos.P());
345  vneg.SetE(vneg.P());
346  TLorentzVector invar = vpos + vneg;
347  upsilon_pair->mass = invar.M();
348 
349  // cuts
350  const bool eta_pos_cut = fabs(gvpos.Eta()) < 1.0;
351  const bool eta_neg_cut = fabs(gvneg.Eta()) < 1.0;
352  const bool reco_upsilon_mass = fabs(upsilon_pair->mass - upsilon_mass) < 0.2; // two sigma cuts
353  upsilon_pair->good_upsilon = eta_pos_cut and eta_neg_cut and reco_upsilon_mass;
354 
355  if (not upsilon_pair->good_upsilon)
356  {
357  return Fun4AllReturnCodes::DISCARDEVENT;
358  }
359  }
360 
361  return Fun4AllReturnCodes::EVENT_OK;
362 }
363 
364 void EMCalAna::eval_trk(PHG4Particle *g4particle, EMCalTrk &trk,
365  PHCompositeNode *topNode)
366 {
367  assert(g4particle);
368 
369  if (!_eval_stack)
370  {
371  _eval_stack = new SvtxEvalStack(topNode);
372  _eval_stack->set_strict(false);
373  }
374  // SvtxVertexEval* vertexeval = _eval_stack->get_vertex_eval();
375  SvtxTrackEval *trackeval = _eval_stack->get_track_eval();
376  // SvtxClusterEval* clustereval = _eval_stack->get_cluster_eval();
377  // SvtxHitEval* hiteval = _eval_stack->get_hit_eval();
378  SvtxTruthEval *trutheval = _eval_stack->get_truth_eval();
379 
380  PHG4TruthInfoContainer *truthinfo =
381  findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
382  assert(truthinfo);
383 
384  SvtxClusterMap *clustermap = findNode::getClass<SvtxClusterMap>(topNode,
385  "SvtxClusterMap");
386  assert(clustermap);
387 
388  int gtrackID = g4particle->get_track_id();
389  int gflavor = g4particle->get_pid();
390 
391  std::set<PHG4Hit *> g4hits = trutheval->all_truth_hits(g4particle);
392  int ng4hits = g4hits.size();
393  float gpx = g4particle->get_px();
394  float gpy = g4particle->get_py();
395  float gpz = g4particle->get_pz();
396 
397  PHG4VtxPoint *vtx = trutheval->get_vertex(g4particle);
398  float gvx = vtx->get_x();
399  float gvy = vtx->get_y();
400  float gvz = vtx->get_z();
401 
402  float gfpx = NAN;
403  float gfpy = NAN;
404  float gfpz = NAN;
405  float gfx = NAN;
406  float gfy = NAN;
407  float gfz = NAN;
408 
409  PHG4Hit *outerhit = trutheval->get_outermost_truth_hit(g4particle);
410  if (outerhit)
411  {
412  gfpx = outerhit->get_px(1);
413  gfpy = outerhit->get_py(1);
414  gfpz = outerhit->get_pz(1);
415  gfx = outerhit->get_x(1);
416  gfy = outerhit->get_y(1);
417  gfz = outerhit->get_z(1);
418  }
419 
420  int gembed = trutheval->get_embed(g4particle);
421 
422  SvtxTrack *track = trackeval->best_track_from(g4particle);
423 
424  float trackID = NAN;
425  float charge = NAN;
426  float quality = NAN;
427  float chisq = NAN;
428  float ndf = NAN;
429  float nhits = NAN;
430  unsigned int layers = 0x0;
431  float dca = NAN;
432  float dca2d = NAN;
433  float dca2dsigma = NAN;
434  float px = NAN;
435  float py = NAN;
436  float pz = NAN;
437  float pcax = NAN;
438  float pcay = NAN;
439  float pcaz = NAN;
440 
441  int nfromtruth = -1;
442 
443  if (track)
444  {
445  trackID = track->get_id();
446  charge = track->get_charge();
447  quality = track->get_quality();
448  chisq = track->get_chisq();
449  ndf = track->get_ndf();
450  nhits = track->size_clusters();
451 
452  for (SvtxTrack::ConstClusterIter iter = track->begin_clusters();
453  iter != track->end_clusters(); ++iter)
454  {
455  unsigned int cluster_id = *iter;
456  SvtxCluster *cluster = clustermap->get(cluster_id);
457  unsigned int layer = cluster->get_layer();
458  if (layer < 32)
459  layers |= (0x1 << layer);
460  }
461 
462  dca = track->get_dca();
463  dca2d = track->get_dca2d();
464  dca2dsigma = track->get_dca2d_error();
465  px = track->get_px();
466  py = track->get_py();
467  pz = track->get_pz();
468  pcax = track->get_x();
469  pcay = track->get_y();
470  pcaz = track->get_z();
471 
472  nfromtruth = trackeval->get_nclusters_contribution(track, g4particle);
473 
474  trk.presdphi = track->get_cal_dphi(SvtxTrack::PRES);
475  trk.presdeta = track->get_cal_deta(SvtxTrack::PRES);
476  trk.prese3x3 = track->get_cal_energy_3x3(SvtxTrack::PRES);
477  trk.prese = track->get_cal_cluster_e(SvtxTrack::PRES);
478 
479  trk.cemcdphi = track->get_cal_dphi(SvtxTrack::CEMC);
480  trk.cemcdeta = track->get_cal_deta(SvtxTrack::CEMC);
481  trk.cemce3x3 = track->get_cal_energy_3x3(SvtxTrack::CEMC);
482  trk.cemce = track->get_cal_cluster_e(SvtxTrack::CEMC);
483 
484  trk.hcalindphi = track->get_cal_dphi(SvtxTrack::HCALIN);
485  trk.hcalindeta = track->get_cal_deta(SvtxTrack::HCALIN);
486  trk.hcaline3x3 = track->get_cal_energy_3x3(SvtxTrack::HCALIN);
487  trk.hcaline = track->get_cal_cluster_e(SvtxTrack::HCALIN);
488 
489  trk.hcaloutdphi = track->get_cal_dphi(SvtxTrack::HCALOUT);
490  trk.hcaloutdeta = track->get_cal_deta(SvtxTrack::HCALOUT);
491  trk.hcaloute3x3 = track->get_cal_energy_3x3(SvtxTrack::HCALOUT);
492  trk.hcaloute = track->get_cal_cluster_e(SvtxTrack::HCALOUT);
493 
494  eval_trk_proj( //
495  /*SvtxTrack **/ track, //
496  /*EMCalTrk &*/ trk,
497  /*enu_calo*/ kCEMC, //
498  /*const double */ gvz,
499  /*PHCompositeNode **/ topNode //
500  );
501  eval_trk_proj( //
502  /*SvtxTrack **/ track, //
503  /*EMCalTrk &*/ trk,
504  /*enu_calo*/ kHCALIN, //
505  /*const double */ gvz,
506  /*PHCompositeNode **/ topNode //
507  );
508  }
509 
510  trk.gtrackID = gtrackID;
511  trk.gflavor = gflavor;
512  trk.ng4hits = ng4hits;
513  trk.gpx = gpx;
514  trk.gpy = gpy;
515  trk.gpz = gpz;
516  trk.gvx = gvx;
517  trk.gvy = gvy;
518  trk.gvz = gvz;
519  trk.gfpx = gfpx;
520  trk.gfpy = gfpy;
521  trk.gfpz = gfpz;
522  trk.gfx = gfx;
523  trk.gfy = gfy;
524  trk.gfz = gfz;
525  trk.gembed = gembed;
526  // trk.gprimary = gprimary;
527  trk.trackID = trackID;
528  trk.px = px;
529  trk.py = py;
530  trk.pz = pz;
531  trk.charge = charge;
532  trk.quality = quality;
533  trk.chisq = chisq;
534  trk.ndf = ndf;
535  trk.nhits = nhits;
536  trk.layers = layers;
537  trk.dca = dca;
538  trk.dca2d = dca2d;
539  trk.dca2dsigma = dca2dsigma;
540  trk.pcax = pcax;
541  trk.pcay = pcay;
542  trk.pcaz = pcaz;
543  trk.nfromtruth = nfromtruth;
544 }
545 
546 void EMCalAna::eval_trk_proj(
547  //
548  SvtxTrack *track, //
549  EMCalTrk &trk, EMCalAna::enu_calo calo_id, const double gvz,
550  PHCompositeNode *topNode //
551  )
552 // Track projections
553 {
554  string detector = "InvalidDetector";
555  double radius = 110;
556  if (calo_id == kCEMC)
557  {
558  detector = "CEMC";
559  radius = 105;
560  }
561  if (calo_id == kHCALIN)
562  {
563  detector = "HCALIN";
564  radius = 125;
565  }
566 
567  if (!_eval_stack)
568  {
569  _eval_stack = new SvtxEvalStack(topNode);
570  _eval_stack->set_strict(false);
571  }
572 
573  // pull the tower geometry
574  string towergeonodename = "TOWERGEOM_" + detector;
575  RawTowerGeomContainer *cemc_towergeo = findNode::getClass<
576  RawTowerGeomContainer>(topNode, towergeonodename.c_str());
577  assert(cemc_towergeo);
578 
579  if (verbosity > 1)
580  {
581  cemc_towergeo->identify();
582  }
583  // pull the towers
584  string towernodename = "TOWER_CALIB_" + detector;
585  RawTowerContainer *cemc_towerList = findNode::getClass<RawTowerContainer>(
586  topNode, towernodename.c_str());
587  assert(cemc_towerList);
588 
589  if (verbosity > 3)
590  {
591  cout << __PRETTY_FUNCTION__ << " - info - handling track track p = ("
592  << track->get_px() << ", " << track->get_py() << ", "
593  << track->get_pz() << "), v = (" << track->get_x() << ", "
594  << track->get_z() << ", " << track->get_z() << ")"
595  << " size_states "
596  << track->size_states() << endl;
597 
598  // for (SvtxTrack::ConstStateIter iter = track->begin_states();
599  // iter != track->end_states(); ++iter)
600  // {
601  // const SvtxTrackState* state = iter->second;
602  // double hitx = state->get_x();
603  // double hity = state->get_y();
604  //
605  // cout << __PRETTY_FUNCTION__ << " - info - track projection : v="
606  // << hitx << ", " << hity <<" p = "<< state->get_px()<<", "<< state->get_py() << endl;
607  // }
608  // track->empty_states();
609 
610  {
611  std::vector<double> point;
612  point.assign(3, -9999.);
613  _hough.projectToRadius(track, _magfield, 10, point);
614 
615  double x = point[0];
616  double y = point[1];
617  double z = point[2];
618  double phi = atan2(y, x);
619  double eta = asinh(z / sqrt(x * x + y * y));
620  cout << __PRETTY_FUNCTION__ << " - info - handling track proj. (" << x
621  << ", " << y << ", " << z << ")"
622  << ", eta " << eta << ", phi"
623  << phi << endl;
624  }
625  {
626  std::vector<double> point;
627  point.assign(3, -9999.);
628  _hough.projectToRadius(track, _magfield, 20, point);
629  double x = point[0];
630  double y = point[1];
631  double z = point[2];
632  double phi = atan2(y, x);
633  double eta = asinh(z / sqrt(x * x + y * y));
634  cout << __PRETTY_FUNCTION__ << " - info - handling track proj. (" << x
635  << ", " << y << ", " << z << ")"
636  << ", eta " << eta << ", phi"
637  << phi << endl;
638  }
639  {
640  std::vector<double> point;
641  point.assign(3, -9999.);
642  _hough.projectToRadius(track, _magfield, 30, point);
643  double x = point[0];
644  double y = point[1];
645  double z = point[2];
646  double phi = atan2(y, x);
647  double eta = asinh(z / sqrt(x * x + y * y));
648  cout << __PRETTY_FUNCTION__ << " - info - handling track proj. (" << x
649  << ", " << y << ", " << z << ")"
650  << ", eta " << eta << ", phi"
651  << phi << endl;
652  }
653  {
654  std::vector<double> point;
655  point.assign(3, -9999.);
656  _hough.projectToRadius(track, _magfield, 40, point);
657  double x = point[0];
658  double y = point[1];
659  double z = point[2];
660  double phi = atan2(y, x);
661  double eta = asinh(z / sqrt(x * x + y * y));
662  cout << __PRETTY_FUNCTION__ << " - info - handling track proj. (" << x
663  << ", " << y << ", " << z << ")"
664  << ", eta " << eta << ", phi"
665  << phi << endl;
666  }
667  {
668  std::vector<double> point;
669  point.assign(3, -9999.);
670  _hough.projectToRadius(track, _magfield, 50, point);
671  double x = point[0];
672  double y = point[1];
673  double z = point[2];
674  double phi = atan2(y, x);
675  double eta = asinh(z / sqrt(x * x + y * y));
676  cout << __PRETTY_FUNCTION__ << " - info - handling track proj. (" << x
677  << ", " << y << ", " << z << ")"
678  << ", eta " << eta << ", phi"
679  << phi << endl;
680  }
681  {
682  std::vector<double> point;
683  point.assign(3, -9999.);
684  _hough.projectToRadius(track, _magfield, 60, point);
685  double x = point[0];
686  double y = point[1];
687  double z = point[2];
688  double phi = atan2(y, x);
689  double eta = asinh(z / sqrt(x * x + y * y));
690  cout << __PRETTY_FUNCTION__ << " - info - handling track proj. (" << x
691  << ", " << y << ", " << z << ")"
692  << ", eta " << eta << ", phi"
693  << phi << endl;
694  }
695  {
696  std::vector<double> point;
697  point.assign(3, -9999.);
698  _hough.projectToRadius(track, _magfield, 70, point);
699  double x = point[0];
700  double y = point[1];
701  double z = point[2];
702  double phi = atan2(y, x);
703  double eta = asinh(z / sqrt(x * x + y * y));
704  cout << __PRETTY_FUNCTION__ << " - info - handling track proj. (" << x
705  << ", " << y << ", " << z << ")"
706  << ", eta " << eta << ", phi"
707  << phi << endl;
708  }
709  {
710  std::vector<double> point;
711  point.assign(3, -9999.);
712  _hough.projectToRadius(track, _magfield, 80, point);
713  double x = point[0];
714  double y = point[1];
715  double z = point[2];
716  double phi = atan2(y, x);
717  double eta = asinh(z / sqrt(x * x + y * y));
718  cout << __PRETTY_FUNCTION__ << " - info - handling track proj. (" << x
719  << ", " << y << ", " << z << ")"
720  << ", eta " << eta << ", phi"
721  << phi << endl;
722  }
723  {
724  std::vector<double> point;
725  point.assign(3, -9999.);
726  _hough.projectToRadius(track, _magfield, 100, point);
727  double x = point[0];
728  double y = point[1];
729  double z = point[2];
730  double phi = atan2(y, x);
731  double eta = asinh(z / sqrt(x * x + y * y));
732  cout << __PRETTY_FUNCTION__ << " - info - handling track proj. (" << x
733  << ", " << y << ", " << z << ")"
734  << ", eta " << eta << ", phi"
735  << phi << endl;
736  }
737  {
738  std::vector<double> point;
739  point.assign(3, -9999.);
740  _hough.projectToRadius(track, _magfield, 1777, point);
741  double x = point[0];
742  double y = point[1];
743  double z = point[2];
744  double phi = atan2(y, x);
745  double eta = asinh(z / sqrt(x * x + y * y));
746  cout << __PRETTY_FUNCTION__ << " - info - handling track proj. (" << x
747  << ", " << y << ", " << z << ")"
748  << ", eta " << eta << ", phi"
749  << phi << endl;
750  }
751  }
752 
753  // curved tracks inside mag field
754  // straight projections thereafter
755  std::vector<double> point;
756  point.assign(3, -9999.);
757  _hough.projectToRadius(track, _magfield, radius, point);
758 
759  if (std::isnan(point[0]) or std::isnan(point[1]) or std::isnan(point[2]))
760  {
761  cout << __PRETTY_FUNCTION__ << "::" << Name()
762  << " - Error - track extrapolation failure:";
763  track->identify();
764  return;
765  }
766 
767  assert(not std::isnan(point[0]));
768  assert(not std::isnan(point[1]));
769  assert(not std::isnan(point[2]));
770 
771  double x = point[0];
772  double y = point[1];
773  // double z = point[2] + gvz - track->get_z();
774  double z = point[2];
775 
776  double phi = atan2(y, x);
777  double eta = asinh(z / sqrt(x * x + y * y));
778 
779  // calculate 3x3 tower energy
780  int binphi = cemc_towergeo->get_phibin(phi);
781  int bineta = cemc_towergeo->get_etabin(eta);
782 
783  double etabin_width = cemc_towergeo->get_etabounds(bineta).second - cemc_towergeo->get_etabounds(bineta).first;
784  if (bineta > 1 and bineta < cemc_towergeo->get_etabins() - 1)
785  etabin_width = (cemc_towergeo->get_etacenter(bineta + 1) - cemc_towergeo->get_etacenter(bineta - 1)) / 2.;
786 
787  double phibin_width = cemc_towergeo->get_phibounds(binphi).second - cemc_towergeo->get_phibounds(binphi).first;
788 
789  assert(etabin_width > 0);
790  assert(phibin_width > 0);
791 
792  const double etabin_shift = (cemc_towergeo->get_etacenter(bineta) - eta) / etabin_width;
793  const double phibin_shift = (fmod(
794  cemc_towergeo->get_phicenter(binphi) - phi + 5 * M_PI, 2 * M_PI) -
795  M_PI) /
796  phibin_width;
797 
798  if (verbosity > 1)
799  cout << __PRETTY_FUNCTION__ << " - info - handling track proj. (" << x
800  << ", " << y << ", " << z << ")"
801  << ", eta " << eta << ", phi" << phi
802  << ", bineta " << bineta << " - ["
803  << cemc_towergeo->get_etabounds(bineta).first << ", "
804  << cemc_towergeo->get_etabounds(bineta).second << "], etabin_shift = "
805  << etabin_shift << ", binphi " << binphi << " - ["
806  << cemc_towergeo->get_phibounds(binphi).first << ", "
807  << cemc_towergeo->get_phibounds(binphi).second << "], phibin_shift = "
808  << phibin_shift << endl;
809 
810  const int bin_search_range = (trk.Max_N_Tower - 1) / 2;
811  for (int iphi = binphi - bin_search_range; iphi <= binphi + bin_search_range;
812  ++iphi)
813 
814  for (int ieta = bineta - bin_search_range;
815  ieta <= bineta + bin_search_range; ++ieta)
816  {
817  // wrap around
818  int wrapphi = iphi;
819  if (wrapphi < 0)
820  {
821  wrapphi = cemc_towergeo->get_phibins() + wrapphi;
822  }
823  if (wrapphi >= cemc_towergeo->get_phibins())
824  {
825  wrapphi = wrapphi - cemc_towergeo->get_phibins();
826  }
827 
828  // edges
829  if (ieta < 0)
830  continue;
831  if (ieta >= cemc_towergeo->get_etabins())
832  continue;
833 
834  if (verbosity > 1)
835  cout << __PRETTY_FUNCTION__ << " - info - processing tower"
836  << ", bineta " << ieta << " - ["
837  << cemc_towergeo->get_etabounds(ieta).first << ", "
838  << cemc_towergeo->get_etabounds(ieta).second << "]"
839  << ", binphi "
840  << wrapphi << " - ["
841  << cemc_towergeo->get_phibounds(wrapphi).first << ", "
842  << cemc_towergeo->get_phibounds(wrapphi).second << "]" << endl;
843 
844  RawTower *tower = cemc_towerList->getTower(ieta, wrapphi);
845  double energy = 0;
846  if (tower)
847  {
848  if (verbosity > 1)
849  cout << __PRETTY_FUNCTION__ << " - info - tower " << ieta << " "
850  << wrapphi << " energy = " << tower->get_energy() << endl;
851 
852  energy = tower->get_energy();
853  }
854 
855  if (calo_id == kCEMC)
856  {
857  trk.cemc_ieta //
858  [ieta - (bineta - bin_search_range)] //
859  [iphi - (binphi - bin_search_range)] //
860  = ieta - bineta + etabin_shift;
861  trk.cemc_iphi //
862  [ieta - (bineta - bin_search_range)] //
863  [iphi - (binphi - bin_search_range)] //
864  = iphi - binphi + phibin_shift;
865  trk.cemc_energy //
866  [ieta - (bineta - bin_search_range)] //
867  [iphi - (binphi - bin_search_range)] //
868  = energy;
869  }
870  if (calo_id == kHCALIN)
871  {
872  trk.hcalin_ieta //
873  [ieta - (bineta - bin_search_range)] //
874  [iphi - (binphi - bin_search_range)] //
875  = ieta - bineta + etabin_shift;
876  trk.hcalin_iphi //
877  [ieta - (bineta - bin_search_range)] //
878  [iphi - (binphi - bin_search_range)] //
879  = iphi - binphi + phibin_shift;
880  trk.hcalin_energy //
881  [ieta - (bineta - bin_search_range)] //
882  [iphi - (binphi - bin_search_range)] //
883  = energy;
884  }
885 
886  } // for (int ieta = bineta-1; ieta < bineta+2; ++ieta) {
887 
888 } // // Track projections
889 
890 int EMCalAna::Init_Trk(PHCompositeNode *topNode)
891 {
892  // static const int BUFFER_SIZE = 32000;
893  // _T_EMCalTrk = new TTree("T_EMCalTrk", "T_EMCalTrk");
894  //
895  // _T_EMCalTrk->Branch("trk.", _trk.ClassName(), &_trk, BUFFER_SIZE);
896  // _T_EMCalTrk->Branch("event", &_ievent, "event/I", BUFFER_SIZE);
897 
898  return Fun4AllReturnCodes::EVENT_OK;
899 }
900 
901 int EMCalAna::process_event_Trk(PHCompositeNode *topNode)
902 {
903  if (verbosity > 2)
904  cout << "EMCalAna::process_event_Trk() entered" << endl;
905 
906  if (!_eval_stack)
907  {
908  _eval_stack = new SvtxEvalStack(topNode);
909  _eval_stack->set_strict(false);
910  }
911 
912  _trk = findNode::getClass<EMCalTrk>(topNode, "EMCalTrk");
913  assert(_trk);
914 
915  PHG4TruthInfoContainer *truthinfo =
916  findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
917  assert(truthinfo);
918  PHG4TruthInfoContainer::ConstRange range =
919  truthinfo->GetPrimaryParticleRange();
920 
921  assert(range.first != range.second);
922  // make sure there is at least one primary partcle
923 
924  eval_trk(range.first->second, *_trk, topNode);
925 
926  return Fun4AllReturnCodes::EVENT_OK;
927 }
928 
929 Fun4AllHistoManager *
930 EMCalAna::get_HistoManager()
931 {
932  Fun4AllServer *se = Fun4AllServer::instance();
933  Fun4AllHistoManager *hm = se->getHistoManager("EMCalAna_HISTOS");
934 
935  if (not hm)
936  {
937  cout
938  << "EMCalAna::get_HistoManager - Making Fun4AllHistoManager EMCalAna_HISTOS"
939  << endl;
940  hm = new Fun4AllHistoManager("EMCalAna_HISTOS");
941  se->registerHistoManager(hm);
942  }
943 
944  assert(hm);
945 
946  return hm;
947 }
948 
949 int EMCalAna::Init_SF(PHCompositeNode *topNode)
950 {
951  Fun4AllHistoManager *hm = get_HistoManager();
952  assert(hm);
953 
954  hm->registerHisto(new TH1F("EMCalAna_h_CEMC_SF", //
955  "_h_CEMC_SF", 1000, 0, .2));
956  hm->registerHisto(new TH1F("EMCalAna_h_HCALIN_SF", //
957  "_h_HCALIN_SF", 1000, 0, .2));
958  hm->registerHisto(new TH1F("EMCalAna_h_HCALOUT_SF", //
959  "_h_HCALOUT_SF", 1000, 0, .2));
960  hm->registerHisto(new TH1F("EMCalAna_h_CEMC_VSF", //
961  "_h_CEMC_VSF", 1000, 0, .2));
962  hm->registerHisto(new TH1F("EMCalAna_h_HCALIN_VSF", //
963  "_h_HCALIN_VSF", 1000, 0, .2));
964  hm->registerHisto(new TH1F("EMCalAna_h_HCALOUT_VSF", //
965  "_h_HCALOUT_VSF", 1000, 0, .2));
966 
967  hm->registerHisto(new TH2F("EMCalAna_h_CEMC_RZ", //
968  "EMCalAna_h_CEMC_RZ;Z (cm);R (cm)", 1200, -300, 300, 600, -000, 300));
969  hm->registerHisto(new TH2F("EMCalAna_h_HCALIN_RZ", //
970  "EMCalAna_h_HCALIN_RZ;Z (cm);R (cm)", 1200, -300, 300, 600, -000, 300));
971  hm->registerHisto(new TH2F("EMCalAna_h_HCALOUT_RZ", //
972  "EMCalAna_h_HCALOUT_RZ;Z (cm);R (cm)", 1200, -300, 300, 600, -000, 300));
973 
974  hm->registerHisto(new TH2F("EMCalAna_h_CEMC_XY", //
975  "EMCalAna_h_CEMC_XY;X (cm);Y (cm)", 1200, -300, 300, 1200, -300, 300));
976  hm->registerHisto(new TH2F("EMCalAna_h_HCALIN_XY", //
977  "EMCalAna_h_HCALIN_XY;X (cm);Y (cm)", 1200, -300, 300, 1200, -300, 300));
978  hm->registerHisto(new TH2F("EMCalAna_h_HCALOUT_XY", //
979  "EMCalAna_h_HCALOUT_XY;X (cm);Y (cm)", 1200, -300, 300, 1200, -300, 300));
980 
981  // block sampling fraction
982  hm->registerHisto(new TH1F("EMCalAna_h_CEMC_BLOCK_EVis", //
983  "EMCalAna_h_CEMC_TOWER_EVis", 100, 0 - .5, 100 - .5));
984  hm->registerHisto(new TH1F("EMCalAna_h_CEMC_BLOCK_ETotal", //
985  "EMCalAna_h_CEMC_BLOCK_ETotal", 100, 0 - .5, 100 - .5));
986 
987  return Fun4AllReturnCodes::EVENT_OK;
988 }
989 
990 int EMCalAna::process_event_SF(PHCompositeNode *topNode)
991 {
992  if (verbosity > 2)
993  cout << "EMCalAna::process_event_SF() entered" << endl;
994 
995  TH1F *h = NULL;
996 
997  Fun4AllHistoManager *hm = get_HistoManager();
998  assert(hm);
999 
1000  double e_hcin = 0.0, e_hcout = 0.0, e_cemc = 0.0;
1001  double ev_hcin = 0.0, ev_hcout = 0.0, ev_cemc = 0.0;
1002  double ea_hcin = 0.0, ea_hcout = 0.0, ea_cemc = 0.0;
1003 
1004  if (_hcalout_hit_container)
1005  {
1006  TH2F *hrz = (TH2F *) hm->getHisto("EMCalAna_h_HCALOUT_RZ");
1007  assert(hrz);
1008  TH2F *hxy = (TH2F *) hm->getHisto("EMCalAna_h_HCALOUT_XY");
1009  assert(hxy);
1010 
1011  PHG4HitContainer::ConstRange hcalout_hit_range =
1012  _hcalout_hit_container->getHits();
1013  for (PHG4HitContainer::ConstIterator hit_iter = hcalout_hit_range.first;
1014  hit_iter != hcalout_hit_range.second; hit_iter++)
1015  {
1016  PHG4Hit *this_hit = hit_iter->second;
1017  assert(this_hit);
1018 
1019  e_hcout += this_hit->get_edep();
1020  ev_hcout += this_hit->get_light_yield();
1021 
1022  const TVector3 hit(this_hit->get_avg_x(), this_hit->get_avg_y(),
1023  this_hit->get_avg_z());
1024  hrz->Fill(hit.Z(), hit.Perp(), this_hit->get_light_yield());
1025  hxy->Fill(hit.X(), hit.Y(), this_hit->get_light_yield());
1026  }
1027  }
1028 
1029  if (_hcalin_hit_container)
1030  {
1031  TH2F *hrz = (TH2F *) hm->getHisto("EMCalAna_h_HCALIN_RZ");
1032  assert(hrz);
1033  TH2F *hxy = (TH2F *) hm->getHisto("EMCalAna_h_HCALIN_XY");
1034  assert(hxy);
1035 
1036  PHG4HitContainer::ConstRange hcalin_hit_range =
1037  _hcalin_hit_container->getHits();
1038  for (PHG4HitContainer::ConstIterator hit_iter = hcalin_hit_range.first;
1039  hit_iter != hcalin_hit_range.second; hit_iter++)
1040  {
1041  PHG4Hit *this_hit = hit_iter->second;
1042  assert(this_hit);
1043 
1044  e_hcin += this_hit->get_edep();
1045  ev_hcin += this_hit->get_light_yield();
1046 
1047  const TVector3 hit(this_hit->get_avg_x(), this_hit->get_avg_y(),
1048  this_hit->get_avg_z());
1049  hrz->Fill(hit.Z(), hit.Perp(), this_hit->get_light_yield());
1050  hxy->Fill(hit.X(), hit.Y(), this_hit->get_light_yield());
1051  }
1052  }
1053 
1054  if (_cemc_hit_container)
1055  {
1056  TH2F *hrz = (TH2F *) hm->getHisto("EMCalAna_h_CEMC_RZ");
1057  assert(hrz);
1058  TH2F *hxy = (TH2F *) hm->getHisto("EMCalAna_h_CEMC_XY");
1059  assert(hxy);
1060 
1061  PHG4HitContainer::ConstRange cemc_hit_range =
1062  _cemc_hit_container->getHits();
1063  for (PHG4HitContainer::ConstIterator hit_iter = cemc_hit_range.first;
1064  hit_iter != cemc_hit_range.second; hit_iter++)
1065  {
1066  PHG4Hit *this_hit = hit_iter->second;
1067  assert(this_hit);
1068 
1069  e_cemc += this_hit->get_edep();
1070  ev_cemc += this_hit->get_light_yield();
1071 
1072  const TVector3 hit(this_hit->get_avg_x(), this_hit->get_avg_y(),
1073  this_hit->get_avg_z());
1074  hrz->Fill(hit.Z(), hit.Perp(), this_hit->get_light_yield());
1075  hxy->Fill(hit.X(), hit.Y(), this_hit->get_light_yield());
1076  }
1077  }
1078 
1079  if (_hcalout_abs_hit_container)
1080  {
1081  PHG4HitContainer::ConstRange hcalout_abs_hit_range =
1082  _hcalout_abs_hit_container->getHits();
1083  for (PHG4HitContainer::ConstIterator hit_iter =
1084  hcalout_abs_hit_range.first;
1085  hit_iter != hcalout_abs_hit_range.second;
1086  hit_iter++)
1087  {
1088  PHG4Hit *this_hit = hit_iter->second;
1089  assert(this_hit);
1090 
1091  ea_hcout += this_hit->get_edep();
1092  }
1093  }
1094 
1095  if (_hcalin_abs_hit_container)
1096  {
1097  PHG4HitContainer::ConstRange hcalin_abs_hit_range =
1098  _hcalin_abs_hit_container->getHits();
1099  for (PHG4HitContainer::ConstIterator hit_iter = hcalin_abs_hit_range.first;
1100  hit_iter != hcalin_abs_hit_range.second; hit_iter++)
1101  {
1102  PHG4Hit *this_hit = hit_iter->second;
1103  assert(this_hit);
1104 
1105  ea_hcin += this_hit->get_edep();
1106  }
1107  }
1108 
1109  if (_cemc_abs_hit_container)
1110  {
1111  PHG4HitContainer::ConstRange cemc_abs_hit_range =
1112  _cemc_abs_hit_container->getHits();
1113  for (PHG4HitContainer::ConstIterator hit_iter = cemc_abs_hit_range.first;
1114  hit_iter != cemc_abs_hit_range.second; hit_iter++)
1115  {
1116  PHG4Hit *this_hit = hit_iter->second;
1117  assert(this_hit);
1118 
1119  ea_cemc += this_hit->get_edep();
1120  }
1121  }
1122 
1123  h = (TH1F *) hm->getHisto("EMCalAna_h_CEMC_SF");
1124  assert(h);
1125  h->Fill(e_cemc / (e_cemc + ea_cemc) + 1e-9);
1126  h = (TH1F *) hm->getHisto("EMCalAna_h_CEMC_VSF");
1127  assert(h);
1128  h->Fill(ev_cemc / (e_cemc + ea_cemc) + 1e-9);
1129 
1130  h = (TH1F *) hm->getHisto("EMCalAna_h_HCALOUT_SF");
1131  assert(h);
1132  h->Fill(e_hcout / (e_hcout + ea_hcout) + 1e-9);
1133  h = (TH1F *) hm->getHisto("EMCalAna_h_HCALOUT_VSF");
1134  assert(h);
1135  h->Fill(ev_hcout / (e_hcout + ea_hcout) + 1e-9);
1136 
1137  h = (TH1F *) hm->getHisto("EMCalAna_h_HCALIN_SF");
1138  assert(h);
1139  h->Fill(e_hcin / (e_hcin + ea_hcin) + 1e-9);
1140  h = (TH1F *) hm->getHisto("EMCalAna_h_HCALIN_VSF");
1141  assert(h);
1142  h->Fill(ev_hcin / (e_hcin + ea_hcin) + 1e-9);
1143 
1144  // block sampling fraction
1145  if (_cemc_hit_container and _cemc_abs_hit_container)
1146  {
1147  TH1 *EMCalAna_h_CEMC_BLOCK_EVis = (TH1 *) hm->getHisto("EMCalAna_h_CEMC_BLOCK_EVis");
1148  assert(EMCalAna_h_CEMC_BLOCK_EVis);
1149  TH1 *EMCalAna_h_CEMC_BLOCK_ETotal = (TH1F *) hm->getHisto("EMCalAna_h_CEMC_BLOCK_ETotal");
1150  assert(EMCalAna_h_CEMC_BLOCK_ETotal);
1151 
1152  const string detector("CEMC");
1153  const string geonodename = "CYLINDERGEOM_" + detector;
1154  const string seggeonodename = "CYLINDERCELLGEOM_" + detector;
1155 
1156  PHG4CylinderGeomContainer *layergeo = findNode::getClass<PHG4CylinderGeomContainer>(topNode, geonodename.c_str());
1157  if (!layergeo)
1158  {
1159  cout << "EMCalAna::process_event_SF - Fatal Error - Could not locate sim geometry node "
1160  << geonodename << endl;
1161  exit(1);
1162  }
1163 
1164  const PHG4CylinderGeom *layergeom_raw = layergeo->GetFirstLayerGeom();
1165  assert(layergeom_raw);
1166  // a special implimentation of PHG4CylinderGeom is required here.
1167  const PHG4CylinderGeom_Spacalv3 *layergeom =
1168  dynamic_cast<const PHG4CylinderGeom_Spacalv3 *>(layergeom_raw);
1169  assert(layergeom);
1170 
1171  PHG4CylinderCellGeomContainer *seggeo = findNode::getClass<PHG4CylinderCellGeomContainer>(topNode, seggeonodename.c_str());
1172  if (!seggeo)
1173  {
1174  cout << "EMCalAna::process_event_SF - Fatal Error - could not locate cell geometry node "
1175  << seggeonodename << endl;
1176  exit(1);
1177  }
1178 
1179  PHG4CylinderCellGeom *geo_raw = seggeo->GetFirstLayerCellGeom();
1180  PHG4CylinderCellGeom_Spacalv1 *geo = dynamic_cast<PHG4CylinderCellGeom_Spacalv1 *>(geo_raw);
1181  assert(geo);
1182 
1183  PHG4HitContainer::ConstRange cemc_hit_range =
1184  _cemc_hit_container->getHits();
1185  for (PHG4HitContainer::ConstIterator hit_iter = cemc_hit_range.first;
1186  hit_iter != cemc_hit_range.second; hit_iter++)
1187  {
1188  PHG4Hit *this_hit = hit_iter->second;
1189  assert(this_hit);
1190  int scint_id = this_hit->get_scint_id();
1191  PHG4CylinderGeom_Spacalv3::scint_id_coder decoder(scint_id);
1192 
1193  std::pair<int, int> tower_z_phi_ID = layergeom->get_tower_z_phi_ID(decoder.tower_ID, decoder.sector_ID);
1194  const int &tower_ID_z = tower_z_phi_ID.first;
1195 
1196  if (Verbosity() >= VERBOSITY_MORE)
1197  {
1198  cout << "EMCalAna::process_event_SF - got scint hit from tower scint_id = " << scint_id
1199  << ", decoder.tower_ID = " << decoder.tower_ID
1200  << ", decoder.sector_ID = " << decoder.sector_ID
1201  << " and decoder tower_ID_z = " << tower_ID_z << endl;
1202  }
1203 
1204  const int etabin = geo->get_etabin_block(tower_ID_z); // block eta bin
1205 
1206  assert(etabin >= 0 and etabin <= 100); // rough range check
1207 
1208  const double e_cemc = this_hit->get_edep();
1209  const double ev_cemc = this_hit->get_light_yield();
1210 
1211  EMCalAna_h_CEMC_BLOCK_EVis->Fill(etabin, ev_cemc);
1212  EMCalAna_h_CEMC_BLOCK_ETotal->Fill(etabin, e_cemc);
1213  }
1214 
1215  PHG4HitContainer::ConstRange cemc_abs_hit_range =
1216  _cemc_abs_hit_container->getHits();
1217  for (PHG4HitContainer::ConstIterator hit_iter = cemc_abs_hit_range.first;
1218  hit_iter != cemc_abs_hit_range.second; hit_iter++)
1219  {
1220  PHG4Hit *this_hit = hit_iter->second;
1221  assert(this_hit);
1222  int scint_id = this_hit->get_scint_id();
1223  PHG4CylinderGeom_Spacalv3::scint_id_coder decoder(scint_id);
1224 
1225  std::pair<int, int> tower_z_phi_ID = layergeom->get_tower_z_phi_ID(decoder.tower_ID, decoder.sector_ID);
1226  const int &tower_ID_z = tower_z_phi_ID.first;
1227 
1228  if (Verbosity() >= VERBOSITY_MORE)
1229  {
1230  cout << "EMCalAna::process_event_SF - got absorber hit from tower scint_id = " << scint_id
1231  << ", decoder.tower_ID = " << decoder.tower_ID
1232  << ", decoder.sector_ID = " << decoder.sector_ID
1233  << " and decoder tower_ID_z = " << tower_ID_z << endl;
1234  }
1235 
1236  if (tower_ID_z > 10)
1237  {
1238  const int etabin = geo->get_etabin_block(tower_ID_z); // block eta bin
1239 
1240  assert(etabin >= 0 and etabin <= 100); // rough range check
1241 
1242  const double ea_cemc = this_hit->get_edep();
1243  EMCalAna_h_CEMC_BLOCK_ETotal->Fill(etabin, ea_cemc);
1244  }
1245  }
1246  }
1247 
1248  return Fun4AllReturnCodes::EVENT_OK;
1249 }
1250 
1251 int EMCalAna::Init_Tower(PHCompositeNode *topNode)
1252 {
1253  Fun4AllHistoManager *hm = get_HistoManager();
1254  assert(hm);
1255 
1256  hm->registerHisto(new TH1F("EMCalAna_h_CEMC_TOWER_1x1", //
1257  "h_CEMC_TOWER_1x1", 5000, 0, 50));
1258  hm->registerHisto(new TH1F("EMCalAna_h_CEMC_TOWER_1x1_max", //
1259  "EMCalAna_h_CEMC_TOWER_1x1_max", 5000, 0, 50));
1260  hm->registerHisto(new TH1F("EMCalAna_h_CEMC_TOWER_1x1_max_trigger_ADC", //
1261  "EMCalAna_h_CEMC_TOWER_1x1_max_trigger_ADC", 5000, 0, 50));
1262 
1263  hm->registerHisto(new TH1F("EMCalAna_h_CEMC_TOWER_2x2", //
1264  "h_CEMC_TOWER_2x2", 5000, 0, 50));
1265  hm->registerHisto(new TH1F("EMCalAna_h_CEMC_TOWER_2x2_max", //
1266  "EMCalAna_h_CEMC_TOWER_2x2_max", 5000, 0, 50));
1267  hm->registerHisto(new TH1F("EMCalAna_h_CEMC_TOWER_2x2_max_trigger_ADC", //
1268  "EMCalAna_h_CEMC_TOWER_2x2_max_trigger_ADC", 5000, 0, 50));
1269  hm->registerHisto(new TH1F("EMCalAna_h_CEMC_TOWER_2x2_slide2_max_trigger_ADC", //
1270  "EMCalAna_h_CEMC_TOWER_2x2_slide2_max_trigger_ADC", 5000, 0, 50));
1271 
1272  hm->registerHisto(new TH1F("EMCalAna_h_CEMC_TOWER_3x3", //
1273  "h_CEMC_TOWER_3x3", 5000, 0, 50));
1274  hm->registerHisto(new TH1F("EMCalAna_h_CEMC_TOWER_3x3_max", //
1275  "EMCalAna_h_CEMC_TOWER_3x3_max", 5000, 0, 50));
1276  hm->registerHisto(new TH1F("EMCalAna_h_CEMC_TOWER_3x3_max_trigger_ADC", //
1277  "EMCalAna_h_CEMC_TOWER_3x3_max_trigger_ADC", 5000, 0, 50));
1278 
1279  hm->registerHisto(new TH1F("EMCalAna_h_CEMC_TOWER_4x4", //
1280  "h_CEMC_TOWER_4x4", 5000, 0, 50));
1281  hm->registerHisto(new TH1F("EMCalAna_h_CEMC_TOWER_4x4_max", //
1282  "EMCalAna_h_CEMC_TOWER_4x4_max", 5000, 0, 50));
1283  hm->registerHisto(new TH1F("EMCalAna_h_CEMC_TOWER_4x4_max_trigger_ADC", //
1284  "EMCalAna_h_CEMC_TOWER_4x4_max_trigger_ADC", 5000, 0, 50));
1285  hm->registerHisto(new TH1F("EMCalAna_h_CEMC_TOWER_4x4_slide2_max_trigger_ADC", //
1286  "EMCalAna_h_CEMC_TOWER_4x4_slide2_max_trigger_ADC", 5000, 0, 50));
1287 
1288  hm->registerHisto(new TH1F("EMCalAna_h_CEMC_TOWER_5x5", //
1289  "h_CEMC_TOWER_4x4", 5000, 0, 50));
1290  hm->registerHisto(new TH1F("EMCalAna_h_CEMC_TOWER_5x5_max", //
1291  "EMCalAna_h_CEMC_TOWER_5x5_max", 5000, 0, 50));
1292  hm->registerHisto(new TH1F("EMCalAna_h_CEMC_TOWER_5x5_max_trigger_ADC", //
1293  "EMCalAna_h_CEMC_TOWER_5x5_max_trigger_ADC", 5000, 0, 50));
1294 
1295  return Fun4AllReturnCodes::EVENT_OK;
1296 }
1297 
1298 int EMCalAna::process_event_Tower(PHCompositeNode *topNode)
1299 {
1300  const string detector("CEMC");
1301 
1302  if (verbosity > 2)
1303  cout << "EMCalAna::process_event_SF() entered" << endl;
1304 
1305  string towernodename = "TOWER_CALIB_" + detector;
1306  // Grab the towers
1307  RawTowerContainer *towers = findNode::getClass<RawTowerContainer>(topNode,
1308  towernodename.c_str());
1309  if (!towers)
1310  {
1311  std::cout << PHWHERE << ": Could not find node " << towernodename.c_str()
1312  << std::endl;
1313  return Fun4AllReturnCodes::DISCARDEVENT;
1314  }
1315  string towergeomnodename = "TOWERGEOM_" + detector;
1316  RawTowerGeomContainer *towergeom = findNode::getClass<RawTowerGeomContainer>(
1317  topNode, towergeomnodename.c_str());
1318  if (!towergeom)
1319  {
1320  cout << PHWHERE << ": Could not find node " << towergeomnodename.c_str()
1321  << endl;
1322  return Fun4AllReturnCodes::ABORTEVENT;
1323  }
1324 
1325  static const double trigger_ADC_bin = 45. / 256.; //8-bit ADC max to 45 GeV
1326  static const int max_size = 5;
1327  map<int, string> size_label;
1328  size_label[1] = "1x1";
1329  size_label[2] = "2x2";
1330  size_label[3] = "3x3";
1331  size_label[4] = "4x4";
1332  size_label[5] = "5x5";
1333  map<int, double> max_energy;
1334  map<int, double> max_energy_trigger_ADC;
1335  map<int, double> slide2_max_energy_trigger_ADC;
1336  map<int, TH1F *> energy_hist_list;
1337  map<int, TH1F *> max_energy_hist_list;
1338  map<int, TH1F *> max_energy_trigger_ADC_hist_list;
1339  map<int, TH1F *> slide2_max_energy_trigger_ADC_hist_list;
1340 
1341  Fun4AllHistoManager *hm = get_HistoManager();
1342  assert(hm);
1343  for (int size = 1; size <= max_size; ++size)
1344  {
1345  max_energy[size] = 0;
1346  max_energy_trigger_ADC[size] = 0;
1347 
1348  TH1F *h = NULL;
1349 
1350  h = (TH1F *) hm->getHisto("EMCalAna_h_CEMC_TOWER_" + size_label[size]);
1351  assert(h);
1352  energy_hist_list[size] = h;
1353  h = (TH1F *) hm->getHisto(
1354  "EMCalAna_h_CEMC_TOWER_" + size_label[size] + "_max");
1355  assert(h);
1356  max_energy_hist_list[size] = h;
1357 
1358  h = (TH1F *) hm->getHisto(
1359  "EMCalAna_h_CEMC_TOWER_" + size_label[size] + "_max_trigger_ADC");
1360  assert(h);
1361  max_energy_trigger_ADC_hist_list[size] = h;
1362 
1363  if (size == 2 or size == 4)
1364  {
1365  // sliding window made from 2x2 sums
1366  slide2_max_energy_trigger_ADC[size] = 0;
1367 
1368  h = (TH1F *) hm->getHisto(
1369  "EMCalAna_h_CEMC_TOWER_" + size_label[size] + "_slide2_max_trigger_ADC");
1370  assert(h);
1371  slide2_max_energy_trigger_ADC_hist_list[size] = h;
1372  }
1373  }
1374 
1375  for (int binphi = 0; binphi < towergeom->get_phibins(); ++binphi)
1376  {
1377  for (int bineta = 0; bineta < towergeom->get_etabins(); ++bineta)
1378  {
1379  for (int size = 1; size <= max_size; ++size)
1380  {
1381  double energy = 0;
1382  double energy_trigger_ADC = 0;
1383  double slide2_energy_trigger_ADC = 0;
1384 
1385  for (int iphi = binphi; iphi < binphi + size; ++iphi)
1386  {
1387  for (int ieta = bineta; ieta < bineta + size; ++ieta)
1388  {
1389  if (ieta > towergeom->get_etabins())
1390  continue;
1391 
1392  // wrap around
1393  int wrapphi = iphi;
1394  assert(wrapphi >= 0);
1395  if (wrapphi >= towergeom->get_phibins())
1396  {
1397  wrapphi = wrapphi - towergeom->get_phibins();
1398  }
1399 
1400  RawTower *tower = towers->getTower(ieta, wrapphi);
1401 
1402  if (tower)
1403  {
1404  const double e_intput = tower->get_energy();
1405 
1406  const double e_trigger_ADC = round(
1407  e_intput / trigger_ADC_bin) *
1408  trigger_ADC_bin;
1409 
1410  energy += e_intput;
1411  energy_trigger_ADC += e_trigger_ADC;
1412 
1413  if ((size == 2 or size == 4) and (binphi % 2 == 0) and (bineta % 2 == 0))
1414  {
1415  // sliding window made from 2x2 sums
1416 
1417  slide2_energy_trigger_ADC += e_trigger_ADC;
1418  }
1419  }
1420  }
1421  }
1422 
1423  energy_hist_list[size]->Fill(energy);
1424 
1425  if (energy > max_energy[size])
1426  max_energy[size] = energy;
1427  if (energy_trigger_ADC > max_energy_trigger_ADC[size])
1428  max_energy_trigger_ADC[size] = energy_trigger_ADC;
1429 
1430  if ((size == 2 or size == 4) and (binphi % 2 == 0) and (bineta % 2 == 0))
1431  {
1432  // sliding window made from 2x2 sums
1433 
1434  if (slide2_energy_trigger_ADC > slide2_max_energy_trigger_ADC[size])
1435  slide2_max_energy_trigger_ADC[size] =
1436  slide2_energy_trigger_ADC;
1437  }
1438 
1439  } // for (int size = 1; size <= 4; ++size)
1440  }
1441  }
1442 
1443  for (int size = 1; size <= max_size; ++size)
1444  {
1445  max_energy_hist_list[size]->Fill(max_energy[size]);
1446  max_energy_trigger_ADC_hist_list[size]->Fill(
1447  max_energy_trigger_ADC[size]);
1448 
1449  if (size == 2 or size == 4)
1450  {
1451  // sliding window made from 2x2 sums
1452  slide2_max_energy_trigger_ADC_hist_list[size]->Fill(
1453  slide2_max_energy_trigger_ADC[size]);
1454  }
1455  }
1456 
1457  return Fun4AllReturnCodes::EVENT_OK;
1458 }