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