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