Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Proto2ShowerCalib.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file Proto2ShowerCalib.C
1 #include "Proto2ShowerCalib.h"
2 #include "TemperatureCorrection.h"
3 
4 #include <prototype2/RawTower_Temperature.h>
5 #include <prototype2/RawTower_Prototype2.h>
6 #include <calobase/RawTowerContainer.h>
7 #include <pdbcalbase/PdbParameterMap.h>
8 #include <phparameter/PHParameters.h>
9 #include <ffaobjects/EventHeader.h>
10 
11 #include <fun4all/SubsysReco.h>
12 #include <fun4all/Fun4AllServer.h>
13 #include <fun4all/PHTFileServer.h>
14 #include <phool/PHCompositeNode.h>
15 #include <fun4all/Fun4AllReturnCodes.h>
16 #include <phool/getClass.h>
17 #include <fun4all/Fun4AllHistoManager.h>
18 
19 #include <phool/PHCompositeNode.h>
20 
21 #include <g4main/PHG4TruthInfoContainer.h>
22 #include <g4main/PHG4Particle.h>
23 #include <g4main/PHG4VtxPoint.h>
24 
25 #include <TNtuple.h>
26 #include <TFile.h>
27 #include <TH1F.h>
28 #include <TH2F.h>
29 #include <TVector3.h>
30 #include <TLorentzVector.h>
31 
32 #include <exception>
33 #include <stdexcept>
34 #include <iostream>
35 #include <sstream>
36 #include <vector>
37 #include <set>
38 #include <algorithm>
39 #include <cassert>
40 #include <cmath>
41 
42 using namespace std;
43 
44 ClassImp(Proto2ShowerCalib::Eval_Cluster);
45 ClassImp(Proto2ShowerCalib::Eval_Run);
46 
47 Proto2ShowerCalib::Proto2ShowerCalib(const std::string &filename) :
48  SubsysReco("Proto2ShowerCalib"), _is_sim(false), _filename(filename), _ievent(
49  0)
50 {
51 
52  verbosity = 1;
53 
54  _eval_run.reset();
55  _eval_3x3_raw.reset();
56  _eval_5x5_raw.reset();
57  _eval_3x3_prod.reset();
58  _eval_5x5_prod.reset();
59  _eval_3x3_temp.reset();
60  _eval_5x5_temp.reset();
61  _eval_3x3_recalib.reset();
62  _eval_5x5_recalib.reset();
63 
64  for (int col = 0; col < n_size; ++col)
65  for (int row = 0; row < n_size; ++row)
66  {
67  _recalib_const[make_pair(col, row)] = 0;
68  }
69 
70 }
71 
72 Proto2ShowerCalib::~Proto2ShowerCalib()
73 {
74 }
75 
76 Fun4AllHistoManager *
77 Proto2ShowerCalib::get_HistoManager()
78 {
79 
80  Fun4AllServer *se = Fun4AllServer::instance();
81  Fun4AllHistoManager *hm = se->getHistoManager("Proto2ShowerCalib_HISTOS");
82 
83  if (not hm)
84  {
85  cout
86  << "Proto2ShowerCalib::get_HistoManager - Making Fun4AllHistoManager Proto2ShowerCalib_HISTOS"
87  << endl;
88  hm = new Fun4AllHistoManager("Proto2ShowerCalib_HISTOS");
89  se->registerHistoManager(hm);
90  }
91 
92  assert(hm);
93 
94  return hm;
95 }
96 
97 int
98 Proto2ShowerCalib::InitRun(PHCompositeNode *topNode)
99 {
100  if (verbosity)
101  cout << "Proto2ShowerCalib::InitRun" << endl;
102 
103  _ievent = 0;
104 
105  PHNodeIterator iter(topNode);
106  PHCompositeNode *dstNode = static_cast<PHCompositeNode*>(iter.findFirst(
107  "PHCompositeNode", "DST"));
108  if (!dstNode)
109  {
110  std::cerr << PHWHERE << "DST Node missing, doing nothing." << std::endl;
111  throw runtime_error(
112  "Failed to find DST node in EmcRawTowerBuilder::CreateNodes");
113  }
114 
115  return Fun4AllReturnCodes::EVENT_OK;
116 }
117 
118 int
119 Proto2ShowerCalib::End(PHCompositeNode *topNode)
120 {
121  cout << "Proto2ShowerCalib::End - write to " << _filename << endl;
122  PHTFileServer::get().cd(_filename);
123 
124  Fun4AllHistoManager *hm = get_HistoManager();
125  assert(hm);
126  for (unsigned int i = 0; i < hm->nHistos(); i++)
127  hm->getHisto(i)->Write();
128 
129 // if (_T_EMCalTrk)
130 // _T_EMCalTrk->Write();
131 
132  fdata.close();
133 
134  return Fun4AllReturnCodes::EVENT_OK;
135 }
136 
137 int
138 Proto2ShowerCalib::Init(PHCompositeNode *topNode)
139 {
140 
141  _ievent = 0;
142 
143  cout << "Proto2ShowerCalib::get_HistoManager - Making PHTFileServer "
144  << _filename << endl;
145  PHTFileServer::get().open(_filename, "RECREATE");
146 
147  fdata.open(_filename + ".dat", std::fstream::out);
148 
149  Fun4AllHistoManager *hm = get_HistoManager();
150  assert(hm);
151 
152  TH2F * hCheck_Cherenkov = new TH2F("hCheck_Cherenkov", "hCheck_Cherenkov",
153  1000, -2000, 2000, 5, .5, 5.5);
154  hCheck_Cherenkov->GetYaxis()->SetBinLabel(1, "C1");
155  hCheck_Cherenkov->GetYaxis()->SetBinLabel(2, "C2 in");
156  hCheck_Cherenkov->GetYaxis()->SetBinLabel(3, "C2 out");
157  hCheck_Cherenkov->GetYaxis()->SetBinLabel(4, "C2 sum");
158  hm->registerHisto(hCheck_Cherenkov);
159 
160  TH1F * hNormalization = new TH1F("hNormalization", "hNormalization", 10, .5,
161  10.5);
162  hCheck_Cherenkov->GetXaxis()->SetBinLabel(1, "ALL");
163  hCheck_Cherenkov->GetXaxis()->SetBinLabel(2, "C2-e");
164  hCheck_Cherenkov->GetXaxis()->SetBinLabel(3, "trigger_veto_pass");
165  hCheck_Cherenkov->GetXaxis()->SetBinLabel(4, "valid_hodo_h");
166  hCheck_Cherenkov->GetXaxis()->SetBinLabel(5, "valid_hodo_v");
167  hCheck_Cherenkov->GetXaxis()->SetBinLabel(6, "good_e");
168  hCheck_Cherenkov->GetXaxis()->SetBinLabel(7, "good_data");
169  hm->registerHisto(hNormalization);
170 
171  hm->registerHisto(new TH1F("hCheck_Veto", "hCheck_Veto", 1000, -500, 500));
172  hm->registerHisto(
173  new TH1F("hCheck_Hodo_H", "hCheck_Hodo_H", 1000, -500, 500));
174  hm->registerHisto(
175  new TH1F("hCheck_Hodo_V", "hCheck_Hodo_V", 1000, -500, 500));
176 
177  hm->registerHisto(new TH1F("hBeam_Mom", "hBeam_Mom", 1200, -120, 120));
178 
179  hm->registerHisto(new TH2F("hEoP", "hEoP", 1000, 0, 1.5, 120, .5, 120.5));
180 
181  hm->registerHisto(new TH1F("hTemperature", "hTemperature", 500, 0, 50));
182 
183  // help index files with TChain
184  TTree * T = new TTree("T", "T");
185  assert(T);
186  hm->registerHisto(T);
187 
188  T->Branch("info", &_eval_run);
189  T->Branch("clus_3x3_raw", &_eval_3x3_raw);
190  T->Branch("clus_5x5_raw", &_eval_5x5_raw);
191  T->Branch("clus_3x3_prod", &_eval_3x3_prod);
192  T->Branch("clus_5x5_prod", &_eval_5x5_prod);
193  T->Branch("clus_3x3_temp", &_eval_3x3_temp);
194  T->Branch("clus_5x5_temp", &_eval_5x5_temp);
195  T->Branch("clus_3x3_recalib", &_eval_3x3_recalib);
196  T->Branch("clus_5x5_recalib", &_eval_5x5_recalib);
197 
198  return Fun4AllReturnCodes::EVENT_OK;
199 }
200 
201 int
202 Proto2ShowerCalib::process_event(PHCompositeNode *topNode)
203 {
204 
205  if (verbosity > 2)
206  cout << "Proto2ShowerCalib::process_event() entered" << endl;
207 
208  // init eval objects
209  _eval_run.reset();
210  _eval_3x3_raw.reset();
211  _eval_5x5_raw.reset();
212  _eval_3x3_prod.reset();
213  _eval_5x5_prod.reset();
214  _eval_3x3_temp.reset();
215  _eval_5x5_temp.reset();
216  _eval_3x3_recalib.reset();
217  _eval_5x5_recalib.reset();
218 
219  Fun4AllHistoManager *hm = get_HistoManager();
220  assert(hm);
221 
222  if (not _is_sim)
223  {
224  PdbParameterMap *info = findNode::getClass<PdbParameterMap>(topNode,
225  "RUN_INFO");
226 
227  assert(info);
228 
229  PHParameters run_info_copy("RunInfo");
230  run_info_copy.FillFrom(info);
231 
232  _eval_run.beam_mom = run_info_copy.get_double_param("beam_MTNRG_GeV");
233 
234  TH1F * hBeam_Mom = dynamic_cast<TH1F *>(hm->getHisto("hBeam_Mom"));
235  assert(hBeam_Mom);
236 
237  hBeam_Mom->Fill(_eval_run.beam_mom);
238  }
239 
240  EventHeader* eventheader = findNode::getClass<EventHeader>(topNode,
241  "EventHeader");
242  if (not _is_sim)
243  {
244  assert(eventheader);
245 
246  _eval_run.run = eventheader->get_RunNumber();
247  if (verbosity > 4)
248  cout << __PRETTY_FUNCTION__ << _eval_run.run << endl;
249 
250  _eval_run.event = eventheader->get_EvtSequence();
251  }
252 
253  if (_is_sim)
254  {
255 
256  PHG4TruthInfoContainer* truthInfoList = findNode::getClass<
257  PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
258 
259  assert(truthInfoList);
260 
261  _eval_run.run = -1;
262 
263  const PHG4Particle * p = truthInfoList->GetPrimaryParticleRange().first->second;
264  assert(p);
265 
266  const PHG4VtxPoint * v = truthInfoList->GetVtx(p->get_vtx_id());
267  assert(v);
268 
269  _eval_run.beam_mom = sqrt(
270  p->get_px() * p->get_px() + p->get_py() * p->get_py()
271  + p->get_pz() * p->get_pz());
272  _eval_run.truth_y = v->get_y();
273  _eval_run.truth_z = v->get_z();
274 
275  }
276 
277  // normalization
278  TH1F * hNormalization = dynamic_cast<TH1F *>(hm->getHisto("hNormalization"));
279  assert(hNormalization);
280 
281  hNormalization->Fill("ALL", 1);
282 
283  RawTowerContainer* TOWER_RAW_CEMC = findNode::getClass<RawTowerContainer>(
284  topNode, _is_sim ? "TOWER_RAW_LG_CEMC" : "TOWER_RAW_CEMC");
285  assert(TOWER_RAW_CEMC);
286  RawTowerContainer* TOWER_CALIB_CEMC = findNode::getClass<RawTowerContainer>(
287  topNode, _is_sim ? "TOWER_CALIB_LG_CEMC" : "TOWER_CALIB_CEMC");
288  assert(TOWER_CALIB_CEMC);
289 
290  // other nodes
291  RawTowerContainer* TOWER_CALIB_TRIGGER_VETO = findNode::getClass<
292  RawTowerContainer>(topNode, "TOWER_CALIB_TRIGGER_VETO");
293  if (not _is_sim)
294  {
295  assert(TOWER_CALIB_TRIGGER_VETO);
296  }
297 
298  RawTowerContainer* TOWER_CALIB_HODO_HORIZONTAL = findNode::getClass<
299  RawTowerContainer>(topNode, "TOWER_CALIB_HODO_HORIZONTAL");
300  if (not _is_sim)
301  {
302  assert(TOWER_CALIB_HODO_HORIZONTAL);
303  }
304  RawTowerContainer* TOWER_CALIB_HODO_VERTICAL = findNode::getClass<
305  RawTowerContainer>(topNode, "TOWER_CALIB_HODO_VERTICAL");
306  if (not _is_sim)
307  {
308  assert(TOWER_CALIB_HODO_VERTICAL);
309  }
310 
311  RawTowerContainer* TOWER_TEMPERATURE_EMCAL = findNode::getClass<
312  RawTowerContainer>(topNode, "TOWER_TEMPERATURE_EMCAL");
313  if (not _is_sim)
314  {
315  assert(TOWER_TEMPERATURE_EMCAL);
316  }
317 
318  RawTowerContainer* TOWER_CALIB_C1 = findNode::getClass<RawTowerContainer>(
319  topNode, "TOWER_CALIB_C1");
320  if (not _is_sim)
321  {
322  assert(TOWER_CALIB_C1);
323  }
324  RawTowerContainer* TOWER_CALIB_C2 = findNode::getClass<RawTowerContainer>(
325  topNode, "TOWER_CALIB_C2");
326  if (not _is_sim)
327  {
328  assert(TOWER_CALIB_C2);
329  }
330 
331  // Cherenkov
332  bool cherekov_e = false;
333  if (not _is_sim)
334  {
335  RawTower * t_c2_in = NULL;
336  RawTower * t_c2_out = NULL;
337 
338  assert(eventheader);
339  if (eventheader->get_RunNumber() >= 2105)
340  {
341  t_c2_in = TOWER_CALIB_C2->getTower(10);
342  t_c2_out = TOWER_CALIB_C2->getTower(11);
343  }
344  else
345  {
346  t_c2_in = TOWER_CALIB_C2->getTower(0);
347  t_c2_out = TOWER_CALIB_C2->getTower(1);
348  }
349  assert(t_c2_in);
350  assert(t_c2_out);
351 
352  const double c2_in = t_c2_in->get_energy();
353  const double c2_out = t_c2_out->get_energy();
354  const double c1 = TOWER_CALIB_C1->getTower(0)->get_energy();
355 
356  _eval_run.C2_sum = c2_in + c2_out;
357  cherekov_e = (_eval_run.C2_sum) > 100;
358  hNormalization->Fill("C2-e", cherekov_e);
359 
360  TH2F * hCheck_Cherenkov = dynamic_cast<TH2F *>(hm->getHisto(
361  "hCheck_Cherenkov"));
362  assert(hCheck_Cherenkov);
363  hCheck_Cherenkov->Fill(c1, "C1", 1);
364  hCheck_Cherenkov->Fill(c2_in, "C2 in", 1);
365  hCheck_Cherenkov->Fill(c2_out, "C2 out", 1);
366  hCheck_Cherenkov->Fill(c2_in + c2_out, "C2 sum", 1);
367  }
368 
369  // veto
370  TH1F * hCheck_Veto = dynamic_cast<TH1F *>(hm->getHisto("hCheck_Veto"));
371  assert(hCheck_Veto);
372  bool trigger_veto_pass = true;
373  if (not _is_sim)
374  {
375  auto range = TOWER_CALIB_TRIGGER_VETO->getTowers();
376  for (auto it = range.first; it != range.second; ++it)
377  {
378  RawTower* tower = it->second;
379  assert(tower);
380 
381  hCheck_Veto->Fill(tower->get_energy());
382 
383  if (abs(tower->get_energy()) > 15)
384  trigger_veto_pass = false;
385  }
386  }
387  hNormalization->Fill("trigger_veto_pass", trigger_veto_pass);
388  _eval_run.trigger_veto_pass = trigger_veto_pass;
389 
390  // hodoscope
391  TH1F * hCheck_Hodo_H = dynamic_cast<TH1F *>(hm->getHisto("hCheck_Hodo_H"));
392  assert(hCheck_Hodo_H);
393  int hodo_h_count = 0;
394  if (not _is_sim)
395  {
396  auto range = TOWER_CALIB_HODO_HORIZONTAL->getTowers();
397  for (auto it = range.first; it != range.second; ++it)
398  {
399  RawTower* tower = it->second;
400  assert(tower);
401 
402  hCheck_Hodo_H->Fill(tower->get_energy());
403 
404  if (abs(tower->get_energy()) > 30)
405  {
406  hodo_h_count++;
407  _eval_run.hodo_h = tower->get_id();
408  }
409  }
410  }
411  const bool valid_hodo_h = hodo_h_count == 1;
412  hNormalization->Fill("valid_hodo_h", valid_hodo_h);
413  _eval_run.valid_hodo_h = valid_hodo_h;
414 
415  TH1F * hCheck_Hodo_V = dynamic_cast<TH1F *>(hm->getHisto("hCheck_Hodo_V"));
416  assert(hCheck_Hodo_V);
417  int hodo_v_count = 0;
418  if (not _is_sim)
419  {
420  auto range = TOWER_CALIB_HODO_VERTICAL->getTowers();
421  for (auto it = range.first; it != range.second; ++it)
422  {
423  RawTower* tower = it->second;
424  assert(tower);
425 
426  hCheck_Hodo_V->Fill(tower->get_energy());
427 
428  if (abs(tower->get_energy()) > 30)
429  {
430  hodo_v_count++;
431  _eval_run.hodo_v = tower->get_id();
432  }
433  }
434  }
435  const bool valid_hodo_v = hodo_v_count == 1;
436  _eval_run.valid_hodo_v = valid_hodo_v;
437  hNormalization->Fill("valid_hodo_v", valid_hodo_v);
438 
439  const bool good_e = (valid_hodo_v and valid_hodo_h and cherekov_e
440  and trigger_veto_pass) and (not _is_sim);
441  hNormalization->Fill("good_e", good_e);
442 
443  // simple clustering
444  pair<int, int> max_3x3 = find_max(TOWER_CALIB_CEMC, 3);
445  pair<int, int> max_5x5 = find_max(TOWER_CALIB_CEMC, 5);
446 
447  _eval_3x3_raw.max_col = max_3x3.first;
448  _eval_3x3_raw.max_row = max_3x3.second;
449  _eval_3x3_prod.max_col = max_3x3.first;
450  _eval_3x3_prod.max_row = max_3x3.second;
451  _eval_3x3_temp.max_col = max_3x3.first;
452  _eval_3x3_temp.max_row = max_3x3.second;
453  _eval_3x3_recalib.max_col = max_3x3.first;
454  _eval_3x3_recalib.max_row = max_3x3.second;
455 
456  _eval_5x5_raw.max_col = max_5x5.first;
457  _eval_5x5_raw.max_row = max_5x5.second;
458  _eval_5x5_prod.max_col = max_5x5.first;
459  _eval_5x5_prod.max_row = max_5x5.second;
460  _eval_5x5_temp.max_col = max_5x5.first;
461  _eval_5x5_temp.max_row = max_5x5.second;
462  _eval_5x5_recalib.max_col = max_5x5.first;
463  _eval_5x5_recalib.max_row = max_5x5.second;
464 
465  // tower
466  bool good_temp = true;
467  double sum_energy_calib = 0;
468  double sum_energy_T = 0;
469  TH1F * hTemperature = dynamic_cast<TH1F *>(hm->getHisto("hTemperature"));
470  assert(hTemperature);
471 
472  stringstream sdata;
473 
474  if (good_e)
475  sdata << abs(_eval_run.beam_mom) << "\t";
476 
477  // tower temperature and recalibration
478  {
479  auto range = TOWER_CALIB_CEMC->getTowers();
480  for (auto it = range.first; it != range.second; ++it)
481  {
482 
483  RawTowerDefs::keytype key = it->first;
484  RawTower* tower = it->second;
485  assert(tower);
486 
487  const int col = tower->get_bineta();
488  const int row = tower->get_binphi();
489 
490  if (col < 0 or col >= 8)
491  continue;
492  if (row < 0 or row >= 8)
493  continue;
494 
495  const double energy_calib = tower->get_energy();
496  sum_energy_calib += energy_calib;
497 
498  RawTower* tower_raw = TOWER_RAW_CEMC->getTower(key);
499  assert(tower_raw);
500 
501  double energy_T = 0;
502  if (not _is_sim)
503  {
504  RawTower_Temperature * temp_t =
505  dynamic_cast<RawTower_Temperature *>(TOWER_TEMPERATURE_EMCAL->getTower(
506  tower->get_row(), tower->get_column())); // note swap of col/row in temperature storage
507  assert(temp_t);
508 
509  const double T = temp_t->get_temperature_from_time(
510  eventheader->get_TimeStamp());
511  hTemperature->Fill(T);
512 
513  if (T < 25 or T > 35)
514  good_temp = false;
515 
516  energy_T = TemperatureCorrection::Apply(energy_calib, T);
517  }
518 
519  // recalibration
520  assert(
521  _recalib_const.find(make_pair(col, row)) != _recalib_const.end());
522  const double energy_recalib = energy_T
523  * _recalib_const[make_pair(col, row)];
524 
525  // energy sums
526  sum_energy_T += energy_T;
527 
528  // calibration file
529 // sdata << tower->get_energy() << "\t";
530  // calibration file - only output 5x5 towers
531  if (col >= max_5x5.first - 2 and col <= max_5x5.first + 2
532  and row >= max_5x5.second - 2 and row <= max_5x5.second + 2)
533  {
534  sdata << tower->get_energy() << "\t";
535  }
536  else
537  {
538  sdata << 0 << "\t";
539  }
540 
541  // cluster 3x3
542  if (col >= max_3x3.first - 1 and col <= max_3x3.first + 1)
543  if (row >= max_3x3.second - 1 and row <= max_3x3.second + 1)
544  {
545  // in cluster
546 
547  _eval_3x3_raw.average_col += abs(tower_raw->get_energy()) * col;
548  _eval_3x3_raw.average_row += abs(tower_raw->get_energy()) * row;
549  _eval_3x3_raw.sum_E += abs(tower_raw->get_energy());
550 
551  _eval_3x3_prod.average_col += energy_calib * col;
552  _eval_3x3_prod.average_row += energy_calib * row;
553  _eval_3x3_prod.sum_E += energy_calib;
554 
555  _eval_3x3_temp.average_col += energy_T * col;
556  _eval_3x3_temp.average_row += energy_T * row;
557  _eval_3x3_temp.sum_E += energy_T;
558 
559  _eval_3x3_recalib.average_col += energy_recalib * col;
560  _eval_3x3_recalib.average_row += energy_recalib * row;
561  _eval_3x3_recalib.sum_E += energy_recalib;
562  }
563 
564  // cluster 5x5
565  if (col >= max_5x5.first - 2 and col <= max_5x5.first + 2)
566  if (row >= max_5x5.second - 2 and row <= max_5x5.second + 2)
567  {
568  // in cluster
569 
570  _eval_5x5_raw.average_col += abs(tower_raw->get_energy()) * col;
571  _eval_5x5_raw.average_row += abs(tower_raw->get_energy()) * row;
572  _eval_5x5_raw.sum_E += abs(tower_raw->get_energy());
573 
574  _eval_5x5_prod.average_col += energy_calib * col;
575  _eval_5x5_prod.average_row += energy_calib * row;
576  _eval_5x5_prod.sum_E += energy_calib;
577 
578  _eval_5x5_temp.average_col += energy_T * col;
579  _eval_5x5_temp.average_row += energy_T * row;
580  _eval_5x5_temp.sum_E += energy_T;
581 
582  _eval_5x5_recalib.average_col += energy_recalib * col;
583  _eval_5x5_recalib.average_row += energy_recalib * row;
584  _eval_5x5_recalib.sum_E += energy_recalib;
585  }
586  }
587  }
588 
589  _eval_3x3_raw.reweight_clus_pol();
590  _eval_5x5_raw.reweight_clus_pol();
591  _eval_3x3_prod.reweight_clus_pol();
592  _eval_5x5_prod.reweight_clus_pol();
593  _eval_3x3_temp.reweight_clus_pol();
594  _eval_5x5_temp.reweight_clus_pol();
595  _eval_3x3_recalib.reweight_clus_pol();
596  _eval_5x5_recalib.reweight_clus_pol();
597 
598  const double EoP = sum_energy_T / abs(_eval_run.beam_mom);
599  hNormalization->Fill("good_temp", good_temp);
600 
601  bool good_data = good_e and good_temp;
602  hNormalization->Fill("good_data", good_data);
603 
604  _eval_run.good_temp = good_temp;
605  _eval_run.good_e = good_e;
606  _eval_run.good_data = good_data;
607  _eval_run.sum_energy_T = sum_energy_T;
608  _eval_run.EoP = EoP;
609 
610  // E/p
611  if (good_data)
612  {
613  if (verbosity >= 3)
614  cout << __PRETTY_FUNCTION__ << " sum_energy_calib = "
615  << sum_energy_calib << " sum_energy_T = " << sum_energy_T
616  << " _eval_run.beam_mom = " << _eval_run.beam_mom << endl;
617 
618  TH2F * hEoP = dynamic_cast<TH2F *>(hm->getHisto("hEoP"));
619  assert(hEoP);
620 
621  hEoP->Fill(EoP, abs(_eval_run.beam_mom));
622  }
623 
624  // calibration file
625  if (good_data and abs(_eval_run.beam_mom) >= 4
626  and abs(_eval_run.beam_mom) <= 8)
627  {
628  assert(fdata.is_open());
629 
630  fdata << sdata.str();
631 
632  fdata << endl;
633  }
634 
635  TTree * T = dynamic_cast<TTree *>(hm->getHisto("T"));
636  assert(T);
637  T->Fill();
638 
639  return Fun4AllReturnCodes::EVENT_OK;
640 }
641 
642 pair<int, int>
643 Proto2ShowerCalib::find_max(RawTowerContainer* towers, int cluster_size)
644 {
645  const int clus_edge_size = (cluster_size - 1) / 2;
646  assert(clus_edge_size >= 0);
647 
648  pair<int, int> max(-1000, -1000);
649  double max_e = 0;
650 
651  for (int col = 0; col < n_size; ++col)
652  for (int row = 0; row < n_size; ++row)
653  {
654  double energy = 0;
655 
656  for (int dcol = col - clus_edge_size; dcol <= col + clus_edge_size;
657  ++dcol)
658  for (int drow = row - clus_edge_size; drow <= row + clus_edge_size;
659  ++drow)
660  {
661  if (dcol < 0 or drow < 0)
662  continue;
663 
664  RawTower * t = towers->getTower(dcol, drow);
665  if (t)
666  energy += t->get_energy();
667  }
668 
669  if (energy > max_e)
670  {
671  max = make_pair(col, row);
672  max_e = energy;
673  }
674  }
675 
676  return max;
677 }
678 
679 int
680 Proto2ShowerCalib::LoadRecalibMap(const std::string & file)
681 {
682  if (verbosity)
683  {
684  cout << __PRETTY_FUNCTION__ << " - input recalibration constant from "
685  << file << endl;
686  }
687 
688  ifstream fcalib(file);
689 
690  assert(fcalib.is_open());
691 
692  string line;
693 
694  while (not fcalib.eof())
695  {
696  getline(fcalib, line);
697 
698  if (verbosity > 10)
699  {
700  cout << __PRETTY_FUNCTION__ << " get line " << line << endl;
701  }
702  istringstream sline(line);
703 
704  int col = -1;
705  int row = -1;
706  double calib = 0;
707 
708  sline >> col >> row >> calib;
709 
710  if (not sline.fail())
711  {
712  if (verbosity)
713  {
714  cout << __PRETTY_FUNCTION__ << " - recalibration constant "
715  << col << "," << row << " = " << calib << endl;
716  }
717 
718  _recalib_const[make_pair(col, row)] = calib;
719  }
720 
721  }
722 
723  return _recalib_const.size();
724 }
725