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