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