Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Proto3ShowerCalib.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file Proto3ShowerCalib.C
1 #include "Proto3ShowerCalib.h"
2 #include "TemperatureCorrection.h"
3 
4 #include <prototype3/RawTower_Temperature.h>
5 #include <prototype3/RawTower_Prototype3.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(Proto3ShowerCalib::Eval_Cluster);
45 ClassImp(Proto3ShowerCalib::Eval_Run);
46 
47 Proto3ShowerCalib::Proto3ShowerCalib(const std::string &filename) :
48  SubsysReco("Proto3ShowerCalib"), _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 Proto3ShowerCalib::~Proto3ShowerCalib()
73 {
74 }
75 
76 Fun4AllHistoManager *
77 Proto3ShowerCalib::get_HistoManager()
78 {
79 
80  Fun4AllServer *se = Fun4AllServer::instance();
81  Fun4AllHistoManager *hm = se->getHistoManager("Proto3ShowerCalib_HISTOS");
82 
83  if (not hm)
84  {
85  cout
86  << "Proto3ShowerCalib::get_HistoManager - Making Fun4AllHistoManager Proto3ShowerCalib_HISTOS"
87  << endl;
88  hm = new Fun4AllHistoManager("Proto3ShowerCalib_HISTOS");
89  se->registerHistoManager(hm);
90  }
91 
92  assert(hm);
93 
94  return hm;
95 }
96 
97 int
98 Proto3ShowerCalib::InitRun(PHCompositeNode *topNode)
99 {
100  if (verbosity)
101  cout << "Proto3ShowerCalib::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 Proto3ShowerCalib::End(PHCompositeNode *topNode)
120 {
121  cout << "Proto3ShowerCalib::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 Proto3ShowerCalib::Init(PHCompositeNode *topNode)
139 {
140 
141  _ievent = 0;
142 
143  cout << "Proto3ShowerCalib::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 Proto3ShowerCalib::process_event(PHCompositeNode *topNode)
203 {
204 
205  if (verbosity > 2)
206  cout << "Proto3ShowerCalib::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  _eval_run.beam_2CH_mm = run_info_copy.get_double_param("beam_2CH_mm");
240  _eval_run.beam_2CV_mm = run_info_copy.get_double_param("beam_2CV_mm");
241 
242  }
243 
244  EventHeader* eventheader = findNode::getClass<EventHeader>(topNode,
245  "EventHeader");
246  if (not _is_sim)
247  {
248  assert(eventheader);
249 
250  _eval_run.run = eventheader->get_RunNumber();
251  if (verbosity > 4)
252  cout << __PRETTY_FUNCTION__ << _eval_run.run << endl;
253 
254  _eval_run.event = eventheader->get_EvtSequence();
255  }
256 
257  if (_is_sim)
258  {
259 
260  PHG4TruthInfoContainer* truthInfoList = findNode::getClass<
261  PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
262 
263  assert(truthInfoList);
264 
265  _eval_run.run = -1;
266 
267  const PHG4Particle * p =
268  truthInfoList->GetPrimaryParticleRange().first->second;
269  assert(p);
270 
271  const PHG4VtxPoint * v = truthInfoList->GetVtx(p->get_vtx_id());
272  assert(v);
273 
274  _eval_run.beam_mom = sqrt(
275  p->get_px() * p->get_px() + p->get_py() * p->get_py()
276  + p->get_pz() * p->get_pz());
277  _eval_run.truth_y = v->get_y();
278  _eval_run.truth_z = v->get_z();
279 
280  }
281 
282  // normalization
283  TH1F * hNormalization = dynamic_cast<TH1F *>(hm->getHisto("hNormalization"));
284  assert(hNormalization);
285 
286  hNormalization->Fill("ALL", 1);
287 
288  RawTowerContainer* TOWER_RAW_CEMC = findNode::getClass<RawTowerContainer>(
289  topNode, _is_sim ? "TOWER_RAW_LG_CEMC" : "TOWER_RAW_CEMC");
290  assert(TOWER_RAW_CEMC);
291  RawTowerContainer* TOWER_CALIB_CEMC = findNode::getClass<RawTowerContainer>(
292  topNode, _is_sim ? "TOWER_CALIB_LG_CEMC" : "TOWER_CALIB_CEMC");
293  assert(TOWER_CALIB_CEMC);
294 
295  // other nodes
296  RawTowerContainer* TOWER_CALIB_TRIGGER_VETO = findNode::getClass<
297  RawTowerContainer>(topNode, "TOWER_CALIB_TRIGGER_VETO");
298  if (not _is_sim)
299  {
300  assert(TOWER_CALIB_TRIGGER_VETO);
301  }
302 
303  RawTowerContainer* TOWER_CALIB_HODO_HORIZONTAL = findNode::getClass<
304  RawTowerContainer>(topNode, "TOWER_CALIB_HODO_HORIZONTAL");
305  if (not _is_sim)
306  {
307  assert(TOWER_CALIB_HODO_HORIZONTAL);
308  }
309  RawTowerContainer* TOWER_CALIB_HODO_VERTICAL = findNode::getClass<
310  RawTowerContainer>(topNode, "TOWER_CALIB_HODO_VERTICAL");
311  if (not _is_sim)
312  {
313  assert(TOWER_CALIB_HODO_VERTICAL);
314  }
315 
316  RawTowerContainer* TOWER_TEMPERATURE_EMCAL = findNode::getClass<
317  RawTowerContainer>(topNode, "TOWER_TEMPERATURE_EMCAL");
318  if (not _is_sim)
319  {
320  assert(TOWER_TEMPERATURE_EMCAL);
321  }
322 
323  RawTowerContainer* TOWER_CALIB_C1 = findNode::getClass<RawTowerContainer>(
324  topNode, "TOWER_CALIB_C1");
325  if (not _is_sim)
326  {
327  assert(TOWER_CALIB_C1);
328  }
329  RawTowerContainer* TOWER_CALIB_C2 = findNode::getClass<RawTowerContainer>(
330  topNode, "TOWER_CALIB_C2");
331  if (not _is_sim)
332  {
333  assert(TOWER_CALIB_C2);
334  }
335 
336  // Cherenkov
337  bool cherekov_e = false;
338  if (not _is_sim)
339  {
340  RawTower * t_c2_in = NULL;
341  RawTower * t_c2_out = NULL;
342 
343  assert(eventheader);
344  if (eventheader->get_RunNumber() >= 2105)
345  {
346  t_c2_in = TOWER_CALIB_C2->getTower(10);
347  t_c2_out = TOWER_CALIB_C2->getTower(11);
348  }
349  else
350  {
351  t_c2_in = TOWER_CALIB_C2->getTower(0);
352  t_c2_out = TOWER_CALIB_C2->getTower(1);
353  }
354  assert(t_c2_in);
355  assert(t_c2_out);
356 
357  const double c2_in = t_c2_in->get_energy();
358  const double c2_out = t_c2_out->get_energy();
359  const double c1 = TOWER_CALIB_C1->getTower(0)->get_energy();
360 
361  _eval_run.C2_sum = c2_in + c2_out;
362  cherekov_e = (_eval_run.C2_sum)
363  > (abs(_eval_run.beam_mom) >= 10 ? 100 : 240);
364  hNormalization->Fill("C2-e", cherekov_e);
365 
366  TH2F * hCheck_Cherenkov = dynamic_cast<TH2F *>(hm->getHisto(
367  "hCheck_Cherenkov"));
368  assert(hCheck_Cherenkov);
369  hCheck_Cherenkov->Fill(c1, "C1", 1);
370  hCheck_Cherenkov->Fill(c2_in, "C2 in", 1);
371  hCheck_Cherenkov->Fill(c2_out, "C2 out", 1);
372  hCheck_Cherenkov->Fill(c2_in + c2_out, "C2 sum", 1);
373  }
374 
375  // veto
376  TH1F * hCheck_Veto = dynamic_cast<TH1F *>(hm->getHisto("hCheck_Veto"));
377  assert(hCheck_Veto);
378  bool trigger_veto_pass = true;
379  if (not _is_sim)
380  {
381  auto range = TOWER_CALIB_TRIGGER_VETO->getTowers();
382  for (auto it = range.first; it != range.second; ++it)
383  {
384  RawTower* tower = it->second;
385  assert(tower);
386 
387  hCheck_Veto->Fill(tower->get_energy());
388 
389  if (abs(tower->get_energy()) > 15)
390  trigger_veto_pass = false;
391  }
392  }
393  hNormalization->Fill("trigger_veto_pass", trigger_veto_pass);
394  _eval_run.trigger_veto_pass = trigger_veto_pass;
395 
396  // hodoscope
397  TH1F * hCheck_Hodo_H = dynamic_cast<TH1F *>(hm->getHisto("hCheck_Hodo_H"));
398  assert(hCheck_Hodo_H);
399  int hodo_h_count = 0;
400  if (not _is_sim)
401  {
402  auto range = TOWER_CALIB_HODO_HORIZONTAL->getTowers();
403  for (auto it = range.first; it != range.second; ++it)
404  {
405  RawTower* tower = it->second;
406  assert(tower);
407 
408  hCheck_Hodo_H->Fill(tower->get_energy());
409 
410  if (abs(tower->get_energy()) > 30)
411  {
412  hodo_h_count++;
413  _eval_run.hodo_h = tower->get_id();
414  }
415  }
416  }
417  const bool valid_hodo_h = hodo_h_count == 1;
418  hNormalization->Fill("valid_hodo_h", valid_hodo_h);
419  _eval_run.valid_hodo_h = valid_hodo_h;
420 
421  TH1F * hCheck_Hodo_V = dynamic_cast<TH1F *>(hm->getHisto("hCheck_Hodo_V"));
422  assert(hCheck_Hodo_V);
423  int hodo_v_count = 0;
424  if (not _is_sim)
425  {
426  auto range = TOWER_CALIB_HODO_VERTICAL->getTowers();
427  for (auto it = range.first; it != range.second; ++it)
428  {
429  RawTower* tower = it->second;
430  assert(tower);
431 
432  hCheck_Hodo_V->Fill(tower->get_energy());
433 
434  if (abs(tower->get_energy()) > 30)
435  {
436  hodo_v_count++;
437  _eval_run.hodo_v = tower->get_id();
438  }
439  }
440  }
441  const bool valid_hodo_v = hodo_v_count == 1;
442  _eval_run.valid_hodo_v = valid_hodo_v;
443  hNormalization->Fill("valid_hodo_v", valid_hodo_v);
444 
445  const bool good_e = (valid_hodo_v and valid_hodo_h and cherekov_e
446  and trigger_veto_pass) and (not _is_sim);
447  hNormalization->Fill("good_e", good_e);
448 
449  // simple clustering
450  pair<int, int> max_3x3 = find_max(TOWER_CALIB_CEMC, 3);
451  pair<int, int> max_5x5 = find_max(TOWER_CALIB_CEMC, 5);
452 
453  _eval_3x3_raw.max_col = max_3x3.first;
454  _eval_3x3_raw.max_row = max_3x3.second;
455  _eval_3x3_prod.max_col = max_3x3.first;
456  _eval_3x3_prod.max_row = max_3x3.second;
457  _eval_3x3_temp.max_col = max_3x3.first;
458  _eval_3x3_temp.max_row = max_3x3.second;
459  _eval_3x3_recalib.max_col = max_3x3.first;
460  _eval_3x3_recalib.max_row = max_3x3.second;
461 
462  _eval_5x5_raw.max_col = max_5x5.first;
463  _eval_5x5_raw.max_row = max_5x5.second;
464  _eval_5x5_prod.max_col = max_5x5.first;
465  _eval_5x5_prod.max_row = max_5x5.second;
466  _eval_5x5_temp.max_col = max_5x5.first;
467  _eval_5x5_temp.max_row = max_5x5.second;
468  _eval_5x5_recalib.max_col = max_5x5.first;
469  _eval_5x5_recalib.max_row = max_5x5.second;
470 
471  // tower
472  bool good_temp = true;
473  double sum_energy_calib = 0;
474  double sum_energy_T = 0;
475  TH1F * hTemperature = dynamic_cast<TH1F *>(hm->getHisto("hTemperature"));
476  assert(hTemperature);
477 
478  stringstream sdata;
479 
480  if (good_e)
481  sdata << abs(_eval_run.beam_mom) << "\t";
482 
483  // tower temperature and recalibration
484  {
485  auto range = TOWER_CALIB_CEMC->getTowers();
486  for (auto it = range.first; it != range.second; ++it)
487  {
488 
489  RawTowerDefs::keytype key = it->first;
490  RawTower* tower = it->second;
491  assert(tower);
492 
493  const int col = tower->get_bineta();
494  const int row = tower->get_binphi();
495 
496  if (col < 0 or col >= 8)
497  continue;
498  if (row < 0 or row >= 8)
499  continue;
500 
501  const double energy_calib = tower->get_energy();
502  sum_energy_calib += energy_calib;
503 
504  RawTower* tower_raw = TOWER_RAW_CEMC->getTower(key);
505  assert(tower_raw);
506 
507  double energy_T = energy_calib;
508 // if (not _is_sim)
509 // {
510 // RawTower_Temperature * temp_t =
511 // dynamic_cast<RawTower_Temperature *>(TOWER_TEMPERATURE_EMCAL->getTower(
512 // tower->get_row(), tower->get_column())); // note swap of col/row in temperature storage
513 // assert(temp_t);
514 //
515 // const double T = temp_t->get_temperature_from_time(
516 // eventheader->get_TimeStamp());
517 // hTemperature->Fill(T);
518 //
519 // if (T < 25 or T > 35)
520 // good_temp = false;
521 //
522 // energy_T = TemperatureCorrection::Apply(energy_calib, T);
523 // }
524 
525  // recalibration
526  assert(
527  _recalib_const.find(make_pair(col, row)) != _recalib_const.end());
528  const double energy_recalib = energy_T
529  * _recalib_const[make_pair(col, row)];
530 
531  // energy sums
532  sum_energy_T += energy_T;
533 
534  // calibration file
535 // sdata << tower->get_energy() << "\t";
536  // calibration file - only output 5x5 towers
537  if (col >= max_5x5.first - 1 and col <= max_5x5.first + 1
538  and row >= max_5x5.second - 1 and row <= max_5x5.second + 1)
539  {
540  sdata << tower->get_energy() << "\t";
541  }
542  else
543  {
544  sdata << 0 << "\t";
545  }
546 
547  // cluster 3x3
548  if (col >= max_3x3.first - 1 and col <= max_3x3.first + 1)
549  if (row >= max_3x3.second - 1 and row <= max_3x3.second + 1)
550  {
551  // in cluster
552 
553  _eval_3x3_raw.average_col += abs(tower_raw->get_energy()) * col;
554  _eval_3x3_raw.average_row += abs(tower_raw->get_energy()) * row;
555  _eval_3x3_raw.sum_E += abs(tower_raw->get_energy());
556 
557  _eval_3x3_prod.average_col += energy_calib * col;
558  _eval_3x3_prod.average_row += energy_calib * row;
559  _eval_3x3_prod.sum_E += energy_calib;
560 
561  _eval_3x3_temp.average_col += energy_T * col;
562  _eval_3x3_temp.average_row += energy_T * row;
563  _eval_3x3_temp.sum_E += energy_T;
564 
565  _eval_3x3_recalib.average_col += energy_recalib * col;
566  _eval_3x3_recalib.average_row += energy_recalib * row;
567  _eval_3x3_recalib.sum_E += energy_recalib;
568  }
569 
570  // cluster 5x5
571  if (col >= max_5x5.first - 2 and col <= max_5x5.first + 2)
572  if (row >= max_5x5.second - 2 and row <= max_5x5.second + 2)
573  {
574  // in cluster
575 
576  _eval_5x5_raw.average_col += abs(tower_raw->get_energy()) * col;
577  _eval_5x5_raw.average_row += abs(tower_raw->get_energy()) * row;
578  _eval_5x5_raw.sum_E += abs(tower_raw->get_energy());
579 
580  _eval_5x5_prod.average_col += energy_calib * col;
581  _eval_5x5_prod.average_row += energy_calib * row;
582  _eval_5x5_prod.sum_E += energy_calib;
583 
584  _eval_5x5_temp.average_col += energy_T * col;
585  _eval_5x5_temp.average_row += energy_T * row;
586  _eval_5x5_temp.sum_E += energy_T;
587 
588  _eval_5x5_recalib.average_col += energy_recalib * col;
589  _eval_5x5_recalib.average_row += energy_recalib * row;
590  _eval_5x5_recalib.sum_E += energy_recalib;
591  }
592  }
593  }
594 
595  _eval_3x3_raw.reweight_clus_pol();
596  _eval_5x5_raw.reweight_clus_pol();
597  _eval_3x3_prod.reweight_clus_pol();
598  _eval_5x5_prod.reweight_clus_pol();
599  _eval_3x3_temp.reweight_clus_pol();
600  _eval_5x5_temp.reweight_clus_pol();
601  _eval_3x3_recalib.reweight_clus_pol();
602  _eval_5x5_recalib.reweight_clus_pol();
603 
604  const double EoP = sum_energy_T / abs(_eval_run.beam_mom);
605  hNormalization->Fill("good_temp", good_temp);
606 
607 // bool good_data = good_e and good_temp;
608  bool good_data = good_e;
609  hNormalization->Fill("good_data", good_data);
610 
611  _eval_run.good_temp = good_temp;
612  _eval_run.good_e = good_e;
613  _eval_run.good_data = good_data;
614  _eval_run.sum_energy_T = sum_energy_T;
615  _eval_run.EoP = EoP;
616 
617  // E/p
618  if (good_data)
619  {
620  if (verbosity >= 3)
621  cout << __PRETTY_FUNCTION__ << " sum_energy_calib = "
622  << sum_energy_calib << " sum_energy_T = " << sum_energy_T
623  << " _eval_run.beam_mom = " << _eval_run.beam_mom << endl;
624 
625  TH2F * hEoP = dynamic_cast<TH2F *>(hm->getHisto("hEoP"));
626  assert(hEoP);
627 
628  hEoP->Fill(EoP, abs(_eval_run.beam_mom));
629  }
630 
631  // calibration file
632  if (good_data and abs(_eval_run.beam_mom) >= 4
633  and abs(_eval_run.beam_mom) <= 8
634  and abs(_eval_3x3_raw.average_col - round(_eval_3x3_raw.average_col))
635  < 0.1
636  and abs(_eval_3x3_raw.average_row - round(_eval_3x3_raw.average_row))
637  < 0.1)
638  {
639  assert(fdata.is_open());
640 
641  fdata << sdata.str();
642 
643  fdata << endl;
644  }
645 
646  TTree * T = dynamic_cast<TTree *>(hm->getHisto("T"));
647  assert(T);
648  T->Fill();
649 
650  return Fun4AllReturnCodes::EVENT_OK;
651 }
652 
653 pair<int, int>
654 Proto3ShowerCalib::find_max(RawTowerContainer* towers, int cluster_size)
655 {
656  const int clus_edge_size = (cluster_size - 1) / 2;
657  assert(clus_edge_size >= 0);
658 
659  pair<int, int> max(-1000, -1000);
660  double max_e = 0;
661 
662  for (int col = 0; col < n_size; ++col)
663  for (int row = 0; row < n_size; ++row)
664  {
665  double energy = 0;
666 
667  for (int dcol = col - clus_edge_size; dcol <= col + clus_edge_size;
668  ++dcol)
669  for (int drow = row - clus_edge_size; drow <= row + clus_edge_size;
670  ++drow)
671  {
672  if (dcol < 0 or drow < 0)
673  continue;
674 
675  RawTower * t = towers->getTower(dcol, drow);
676  if (t)
677  energy += t->get_energy();
678  }
679 
680  if (energy > max_e)
681  {
682  max = make_pair(col, row);
683  max_e = energy;
684  }
685  }
686 
687  return max;
688 }
689 
690 int
691 Proto3ShowerCalib::LoadRecalibMap(const std::string & file)
692 {
693  if (verbosity)
694  {
695  cout << __PRETTY_FUNCTION__ << " - input recalibration constant from "
696  << file << endl;
697  }
698 
699  ifstream fcalib(file);
700 
701  assert(fcalib.is_open());
702 
703  string line;
704 
705  while (not fcalib.eof())
706  {
707  getline(fcalib, line);
708 
709  if (verbosity > 10)
710  {
711  cout << __PRETTY_FUNCTION__ << " get line " << line << endl;
712  }
713  istringstream sline(line);
714 
715  int col = -1;
716  int row = -1;
717  double calib = 0;
718 
719  sline >> col >> row >> calib;
720 
721  if (not sline.fail())
722  {
723  if (verbosity)
724  {
725  cout << __PRETTY_FUNCTION__ << " - recalibration constant "
726  << col << "," << row << " = " << calib << endl;
727  }
728 
729  _recalib_const[make_pair(col, row)] = calib;
730  }
731 
732  }
733 
734  return _recalib_const.size();
735 }
736