HCAL/ShowerCalib/Proto2ShowerCalib.C

0001 #include "Proto2ShowerCalib.h"
0002
0003 #include <prototype2/RawTower_Prototype2.h>
0004 #include <calobase/RawTowerContainer.h>
0005 #include <pdbcalbase/PdbParameterMap.h>
0006 #include <phparameter/PHParameters.h>
0007 #include <ffaobjects/EventHeader.h>
0008
0009 #include <fun4all/SubsysReco.h>
0010 #include <fun4all/Fun4AllServer.h>
0011 #include <fun4all/PHTFileServer.h>
0012 #include <phool/PHCompositeNode.h>
0013 #include <fun4all/Fun4AllReturnCodes.h>
0014 #include <phool/getClass.h>
0015 #include <fun4all/Fun4AllHistoManager.h>
0016
0017 #include <phool/PHCompositeNode.h>
0018
0019 #include <TNtuple.h>
0020 #include <TFile.h>
0021 #include <TH1F.h>
0022 #include <TH2F.h>
0023 #include <TVector3.h>
0024 #include <TLorentzVector.h>
0025
0026 #include <exception>
0027 #include <stdexcept>
0028 #include <iostream>
0029 #include <sstream>
0030 #include <vector>
0031 #include <set>
0032 #include <algorithm>
0033 #include <cassert>
0034 #include <cmath>
0035
0036 using namespace std;
0037
0038 ClassImp(Proto2ShowerCalib::HCAL_shower);
0039 ClassImp(Proto2ShowerCalib::Eval_Run);
0040 ClassImp(Proto2ShowerCalib::Time_Samples);
0041
0042 Proto2ShowerCalib::Proto2ShowerCalib(const std::string &filename) :
0043     SubsysReco("Proto2ShowerCalib"), _filename(filename), _ievent(0)
0044 {
0045
0046   verbosity = 1;
0047   _is_sim = false;
0048   _fill_time_samples = false;
0049   _eval_run.reset();
0050   _shower.reset();
0051   _time_samples.reset();
0052
0053   for (int col = 0; col < n_size_emcal; ++col)
0054     for (int row = 0; row < n_size_emcal; ++row)
0055       {
0056         _emcal_recalib_const[make_pair(col, row)] = 0;
0057       }
0058
0059   for (int col = 0; col < n_size_hcalin; ++col)
0060     for (int row = 0; row < n_size_hcalin; ++row)
0061      {
0062         _hcalin_recalib_const[make_pair(col, row)] = 0;
0063         _hcalout_recalib_const[make_pair(col, row)] = 0;
0064      }
0065
0066 }
0067
0068 void Proto2ShowerCalib::fill_time_samples(const bool b)
0069 {
0070   _fill_time_samples = b;
0071 }
0072
0073
0074
0075 Proto2ShowerCalib::~Proto2ShowerCalib()
0076 {
0077 }
0078
0079 Fun4AllHistoManager *
0080 Proto2ShowerCalib::get_HistoManager()
0081 {
0082
0083   Fun4AllServer *se = Fun4AllServer::instance();
0084   Fun4AllHistoManager *hm = se->getHistoManager("Proto2ShowerCalib_HISTOS");
0085
0086   if (not hm)
0087     {
0088       cout
0089           << "Proto2ShowerCalib::get_HistoManager - Making Fun4AllHistoManager Proto2ShowerCalib_HISTOS"
0090           << endl;
0091       hm = new Fun4AllHistoManager("Proto2ShowerCalib_HISTOS");
0092       se->registerHistoManager(hm);
0093     }
0094
0095   assert(hm);
0096
0097   return hm;
0098 }
0099
0100 int
0101 Proto2ShowerCalib::InitRun(PHCompositeNode *topNode)
0102 {
0103   if (verbosity)
0104     cout << "Proto2ShowerCalib::InitRun" << endl;
0105
0106   _ievent = 0;
0107
0108   PHNodeIterator iter(topNode);
0109   PHCompositeNode *dstNode = static_cast<PHCompositeNode*>(iter.findFirst(
0110       "PHCompositeNode", "DST"));
0111   if (!dstNode)
0112     {
0113       std::cerr << PHWHERE << "DST Node missing, doing nothing." << std::endl;
0114       throw runtime_error(
0115           "Failed to find DST node in EmcRawTowerBuilder::CreateNodes");
0116     }
0117
0118   return Fun4AllReturnCodes::EVENT_OK;
0119 }
0120
0121 int
0122 Proto2ShowerCalib::End(PHCompositeNode *topNode)
0123 {
0124   cout << "Proto2ShowerCalib::End - write to " << _filename << endl;
0125   PHTFileServer::get().cd(_filename);
0126
0127   Fun4AllHistoManager *hm = get_HistoManager();
0128   assert(hm);
0129   for (unsigned int i = 0; i < hm->nHistos(); i++)
0130     hm->getHisto(i)->Write();
0131
0132 //  if (_T_EMCalTrk)
0133 //    _T_EMCalTrk->Write();
0134
0135   fdata.close();
0136
0137   return Fun4AllReturnCodes::EVENT_OK;
0138 }
0139
0140 int
0141 Proto2ShowerCalib::Init(PHCompositeNode *topNode)
0142 {
0143
0144   _ievent = 0;
0145
0146   cout << "Proto2ShowerCalib::get_HistoManager - Making PHTFileServer "
0147       << _filename << endl;
0148   PHTFileServer::get().open(_filename, "RECREATE");
0149
0150   fdata.open(_filename + ".dat", std::fstream::out);
0151
0152   Fun4AllHistoManager *hm = get_HistoManager();
0153   assert(hm);
0154
0155   TH2F * hCheck_Cherenkov = new TH2F("hCheck_Cherenkov", "hCheck_Cherenkov",
0156       1000, -2000, 2000, 5, .5, 5.5);
0157   hCheck_Cherenkov->GetYaxis()->SetBinLabel(1, "C1");
0158   hCheck_Cherenkov->GetYaxis()->SetBinLabel(2, "C2 in");
0159   hCheck_Cherenkov->GetYaxis()->SetBinLabel(3, "C2 out");
0160   hCheck_Cherenkov->GetYaxis()->SetBinLabel(4, "C2 sum");
0161   hm->registerHisto(hCheck_Cherenkov);
0162
0163   TH1F * hNormalization = new TH1F("hNormalization", "hNormalization", 10, .5,
0164       10.5);
0165   hCheck_Cherenkov->GetXaxis()->SetBinLabel(1, "ALL");
0166   hCheck_Cherenkov->GetXaxis()->SetBinLabel(2, "C2-e");
0167   hCheck_Cherenkov->GetXaxis()->SetBinLabel(3, "C2-h");
0168   hCheck_Cherenkov->GetXaxis()->SetBinLabel(4, "trigger_veto_pass");
0169   hCheck_Cherenkov->GetXaxis()->SetBinLabel(5, "valid_hodo_h");
0170   hCheck_Cherenkov->GetXaxis()->SetBinLabel(6, "valid_hodo_v");
0171   hCheck_Cherenkov->GetXaxis()->SetBinLabel(7, "good_e");
0172   hCheck_Cherenkov->GetXaxis()->SetBinLabel(8, "good_h");
0173   hm->registerHisto(hNormalization);
0174
0175   hm->registerHisto(new TH1F("hCheck_Veto", "hCheck_Veto", 1000, -500, 500));
0176   hm->registerHisto(
0177       new TH1F("hCheck_Hodo_H", "hCheck_Hodo_H", 1000, -500, 500));
0178   hm->registerHisto(
0179       new TH1F("hCheck_Hodo_V", "hCheck_Hodo_V", 1000, -500, 500));
0180
0181   hm->registerHisto(new TH1F("hBeam_Mom", "hBeam_Mom", 1200, -120, 120));
0182
0183   hm->registerHisto(new TH2F("hEoP", "hEoP", 1000, 0, 1.5, 120, .5, 120.5));
0184
0185   hm->registerHisto(new TH1F("hTemperature", "hTemperature", 500, 0, 50));
0186
0187   // help index files with TChain
0188   TTree * T = new TTree("T", "T");
0189   assert(T);
0190   hm->registerHisto(T);
0191
0192   T->Branch("info", &_eval_run);
0193   T->Branch("shower", &_shower);
0194   if(_fill_time_samples) T->Branch("time", &_time_samples);
0195
0196   return Fun4AllReturnCodes::EVENT_OK;
0197 }
0198
0199 int
0200 Proto2ShowerCalib::process_event(PHCompositeNode *topNode)
0201 {
0202
0203   if (verbosity > 2)
0204     cout << "Proto2ShowerCalib::process_event() entered" << endl;
0205
0206   // init eval objects
0207   _eval_run.reset();
0208   _shower.reset();
0209
0210   Fun4AllHistoManager *hm = get_HistoManager();
0211   assert(hm);
0212
0213   PdbParameterMap *info = findNode::getClass<PdbParameterMap>(topNode,
0214       "RUN_INFO");
0215
0216   //if(!_is_sim) assert(info);
0217
0218   if(info)
0219   {
0220    PHParameters run_info_copy("RunInfo");
0221    run_info_copy.FillFrom(info);
0222
0223    _eval_run.beam_mom = run_info_copy.get_double_param("beam_MTNRG_GeV");
0224
0225    TH1F * hBeam_Mom = dynamic_cast<TH1F *>(hm->getHisto("hBeam_Mom"));
0226    assert(hBeam_Mom);
0227
0228    hBeam_Mom->Fill(_eval_run.beam_mom);
0229   }
0230
0231   EventHeader* eventheader = findNode::getClass<EventHeader>(topNode,
0232       "EventHeader");
0233   //if(!_is_sim) assert(eventheader);
0234
0235   if( eventheader )
0236   {
0237    _eval_run.run = eventheader->get_RunNumber();
0238    if (verbosity > 4)
0239     cout << __PRETTY_FUNCTION__ << _eval_run.run << endl;
0240
0241    _eval_run.event = eventheader->get_EvtSequence();
0242   }
0243   // normalization
0244   TH1F * hNormalization = dynamic_cast<TH1F *>(hm->getHisto("hNormalization"));
0245   assert(hNormalization);
0246
0247   hNormalization->Fill("ALL", 1);
0248
0249   // other nodes
0250   RawTowerContainer* TOWER_CALIB_TRIGGER_VETO = findNode::getClass<
0251       RawTowerContainer>(topNode, "TOWER_CALIB_TRIGGER_VETO");
0252   //if(!_is_sim) assert(TOWER_CALIB_TRIGGER_VETO);
0253
0254   RawTowerContainer* TOWER_CALIB_HODO_HORIZONTAL = findNode::getClass<
0255       RawTowerContainer>(topNode, "TOWER_CALIB_HODO_HORIZONTAL");
0256   //if(!_is_sim) assert(TOWER_CALIB_HODO_HORIZONTAL);
0257
0258   RawTowerContainer* TOWER_CALIB_HODO_VERTICAL = findNode::getClass<
0259       RawTowerContainer>(topNode, "TOWER_CALIB_HODO_VERTICAL");
0260   //if(!_is_sim) assert(TOWER_CALIB_HODO_VERTICAL);
0261
0262  RawTowerContainer* TOWER_RAW_CEMC;
0263   if(!_is_sim)
0264   {
0265    TOWER_RAW_CEMC = findNode::getClass<RawTowerContainer>(
0266       topNode, "TOWER_RAW_CEMC");
0267   }else{
0268    TOWER_RAW_CEMC = findNode::getClass<RawTowerContainer>(
0269       topNode, "TOWER_RAW_LG_CEMC");
0270  }
0271   assert(TOWER_RAW_CEMC);
0272
0273  RawTowerContainer* TOWER_CALIB_CEMC;
0274  if(!_is_sim)
0275  {
0276   TOWER_CALIB_CEMC = findNode::getClass<RawTowerContainer>(
0277       topNode, "TOWER_CALIB_CEMC");
0278  } else {
0279   TOWER_CALIB_CEMC = findNode::getClass<RawTowerContainer>(
0280       topNode, "TOWER_CALIB_LG_CEMC");
0281  }
0282   assert(TOWER_CALIB_CEMC);
0283
0284   RawTowerContainer* TOWER_RAW_LG_HCALIN = 0;
0285   RawTowerContainer* TOWER_RAW_HG_HCALIN = 0;
0286   if(!_is_sim)
0287   {
0288    TOWER_RAW_LG_HCALIN = findNode::getClass<RawTowerContainer>(
0289       topNode, "TOWER_RAW_LG_HCALIN");
0290    TOWER_RAW_HG_HCALIN= findNode::getClass<RawTowerContainer>(
0291       topNode, "TOWER_RAW_HG_HCALIN");
0292    assert(TOWER_RAW_LG_HCALIN);
0293    assert(TOWER_RAW_HG_HCALIN);
0294   }
0295
0296   RawTowerContainer* TOWER_CALIB_HCALIN;
0297   if(!_is_sim)
0298   {
0299    TOWER_CALIB_HCALIN = findNode::getClass<RawTowerContainer>(
0300       topNode, "TOWER_CALIB_LG_HCALIN");
0301   } else {
0302    TOWER_CALIB_HCALIN = findNode::getClass<RawTowerContainer>(
0303       topNode, "TOWER_CALIB_LG_HCALIN");
0304   }
0305   assert(TOWER_CALIB_HCALIN);
0306
0307   RawTowerContainer* TOWER_RAW_LG_HCALOUT = 0;
0308   RawTowerContainer* TOWER_RAW_HG_HCALOUT = 0;
0309   if(!_is_sim)
0310   {
0311    TOWER_RAW_LG_HCALOUT = findNode::getClass<RawTowerContainer>(
0312       topNode, "TOWER_RAW_LG_HCALOUT");
0313    TOWER_RAW_HG_HCALOUT = findNode::getClass<RawTowerContainer>(
0314       topNode, "TOWER_RAW_HG_HCALOUT");
0315    assert(TOWER_RAW_LG_HCALOUT);
0316    assert(TOWER_RAW_HG_HCALOUT);
0317   }
0318
0319   RawTowerContainer* TOWER_CALIB_HCALOUT = findNode::getClass<RawTowerContainer>(
0320       topNode, "TOWER_CALIB_LG_HCALOUT");
0321   assert(TOWER_CALIB_HCALOUT);
0322
0323   RawTowerContainer* TOWER_TEMPERATURE_EMCAL = findNode::getClass<
0324       RawTowerContainer>(topNode, "TOWER_TEMPERATURE_EMCAL");
0325   if(!_is_sim) assert(TOWER_TEMPERATURE_EMCAL);
0326
0327   RawTowerContainer* TOWER_CALIB_C1 = findNode::getClass<RawTowerContainer>(
0328       topNode, "TOWER_CALIB_C1");
0329   //if(!_is_sim) assert(TOWER_CALIB_C1);
0330   RawTowerContainer* TOWER_CALIB_C2 = findNode::getClass<RawTowerContainer>(
0331       topNode, "TOWER_CALIB_C2");
0332   //if(!_is_sim) assert(TOWER_CALIB_C2);
0333
0334   if(!_is_sim && TOWER_CALIB_C2 && TOWER_CALIB_C1 && eventheader)
0335   {
0336   // Cherenkov
0337   RawTower * t_c2_in = NULL;
0338   RawTower * t_c2_out = NULL;
0339
0340   if (eventheader->get_RunNumber() >= 2105)
0341     {
0342       t_c2_in = TOWER_CALIB_C2->getTower(10);
0343       t_c2_out = TOWER_CALIB_C2->getTower(11);
0344     }
0345   else
0346     {
0347       t_c2_in = TOWER_CALIB_C2->getTower(0);
0348       t_c2_out = TOWER_CALIB_C2->getTower(1);
0349     }
0350   assert(t_c2_in);
0351   assert(t_c2_out);
0352
0353   const double c2_in = t_c2_in->get_energy();
0354   const double c2_out = t_c2_out->get_energy();
0355   const double c1 = TOWER_CALIB_C1->getTower(0)->get_energy();
0356
0357   _eval_run.C2_sum = c2_in + c2_out;
0358   _eval_run.C1 = c1;
0359   bool cherekov_e = (_eval_run.C2_sum) > 100;
0360   hNormalization->Fill("C2-e", cherekov_e);
0361
0362   bool cherekov_h = (_eval_run.C2_sum) < 20;
0363   hNormalization->Fill("C2-h", cherekov_h);
0364
0365   TH2F * hCheck_Cherenkov = dynamic_cast<TH2F *>(hm->getHisto(
0366       "hCheck_Cherenkov"));
0367   assert(hCheck_Cherenkov);
0368   hCheck_Cherenkov->Fill(c1, "C1", 1);
0369   hCheck_Cherenkov->Fill(c2_in, "C2 in", 1);
0370   hCheck_Cherenkov->Fill(c2_out, "C2 out", 1);
0371   hCheck_Cherenkov->Fill(c2_in + c2_out, "C2 sum", 1);
0372
0373   // veto
0374   TH1F * hCheck_Veto = dynamic_cast<TH1F *>(hm->getHisto("hCheck_Veto"));
0375   assert(hCheck_Veto);
0376   bool trigger_veto_pass = true;
0377     {
0378       auto range = TOWER_CALIB_TRIGGER_VETO->getTowers();
0379       for (auto it = range.first; it != range.second; ++it)
0380         {
0381           RawTower* tower = it->second;
0382           assert(tower);
0383
0384           hCheck_Veto->Fill(tower->get_energy());
0385
0386           if (abs(tower->get_energy()) > 15)
0387             trigger_veto_pass = false;
0388         }
0389     }
0390   hNormalization->Fill("trigger_veto_pass", trigger_veto_pass);
0391   _eval_run.trigger_veto_pass = trigger_veto_pass;
0392
0393   // hodoscope
0394   TH1F * hCheck_Hodo_H = dynamic_cast<TH1F *>(hm->getHisto("hCheck_Hodo_H"));
0395   assert(hCheck_Hodo_H);
0396   int hodo_h_count = 0;
0397     {
0398       auto range = TOWER_CALIB_HODO_HORIZONTAL->getTowers();
0399       for (auto it = range.first; it != range.second; ++it)
0400         {
0401           RawTower* tower = it->second;
0402           assert(tower);
0403           hCheck_Hodo_H->Fill(tower->get_energy());
0404
0405           if (abs(tower->get_energy()) > 30)
0406             {
0407               hodo_h_count++;
0408               _eval_run.hodo_h = tower->get_id();
0409             }
0410         }
0411     }
0412
0413   const bool valid_hodo_h = hodo_h_count == 1;
0414   hNormalization->Fill("valid_hodo_h", valid_hodo_h);
0415   _eval_run.valid_hodo_h = valid_hodo_h;
0416
0417   TH1F * hCheck_Hodo_V = dynamic_cast<TH1F *>(hm->getHisto("hCheck_Hodo_V"));
0418   assert(hCheck_Hodo_V);
0419   int hodo_v_count = 0;
0420    {
0421       auto range = TOWER_CALIB_HODO_VERTICAL->getTowers();
0422       for (auto it = range.first; it != range.second; ++it)
0423         {
0424           RawTower* tower = it->second;
0425           assert(tower);
0426
0427           hCheck_Hodo_V->Fill(tower->get_energy());
0428
0429           if (abs(tower->get_energy()) > 30)
0430             {
0431               hodo_v_count++;
0432               _eval_run.hodo_v = tower->get_id();
0433             }
0434         }
0435     }
0436   const bool valid_hodo_v = hodo_v_count == 1;
0437   _eval_run.valid_hodo_v = valid_hodo_v;
0438   hNormalization->Fill("valid_hodo_v", valid_hodo_v);
0439
0440   const bool good_e = valid_hodo_v and valid_hodo_h and cherekov_e and trigger_veto_pass;
0441
0442   const bool good_h = valid_hodo_v and valid_hodo_h and cherekov_h and trigger_veto_pass;
0443
0444   hNormalization->Fill("good_e", good_e);
0445   hNormalization->Fill("good_h", good_h);
0446
0447   _eval_run.good_temp = 1;
0448   _eval_run.good_e = good_e;
0449   _eval_run.good_h = good_h;
0450
0451   }
0452   // tower
0453   double emcal_sum_energy_calib = 0;
0454   double emcal_sum_energy_recalib = 0;
0455
0456   double hcalin_sum_energy_calib = 0;
0457   double hcalin_sum_energy_recalib = 0;
0458
0459   double hcalout_sum_energy_calib = 0;
0460   double hcalout_sum_energy_recalib = 0;
0461
0462   stringstream sdata;
0463
0464    //EMCAL towers
0465     {
0466       auto range = TOWER_CALIB_CEMC->getTowers();
0467       for (auto it = range.first; it != range.second; ++it)
0468         {
0469
0470           //RawTowerDefs::keytype key = it->first;
0471           RawTower* tower = it->second;
0472           assert(tower);
0473
0474           const double energy_calib = tower->get_energy();
0475           emcal_sum_energy_calib += energy_calib;
0476
0477           if(_is_sim) continue;
0478           const int col = tower->get_column();
0479           const int row = tower->get_row();
0480
0481           // recalibration
0482           assert(
0483               _emcal_recalib_const.find(make_pair(col, row)) != _emcal_recalib_const.end());
0484           const double energy_recalib = energy_calib
0485               * _emcal_recalib_const[make_pair(col, row)];
0486
0487           // energy sums
0488           emcal_sum_energy_recalib += energy_recalib;
0489          }
0490       }
0491
0492   //HCALIN towers
0493   {
0494       auto range = TOWER_CALIB_HCALIN->getTowers();
0495       for (auto it = range.first; it != range.second; ++it)
0496         {
0497           RawTower* tower = it->second;
0498           assert(tower);
0499
0500           const double energy_calib = tower->get_energy();
0501           hcalin_sum_energy_calib += energy_calib;
0502
0503           if(_is_sim) continue;
0504           const int col = tower->get_column();
0505           const int row = tower->get_row();
0506
0507           assert(
0508               _hcalin_recalib_const.find(make_pair(col, row)) != _hcalin_recalib_const.end());
0509           const double energy_recalib = energy_calib
0510               * _hcalin_recalib_const[make_pair(col, row)];
0511
0512           hcalin_sum_energy_recalib += energy_recalib;
0513          }
0514      }
0515
0516    //HCALOUT towers
0517    {
0518       auto range = TOWER_CALIB_HCALOUT->getTowers();
0519       for (auto it = range.first; it != range.second; ++it)
0520         {
0521           RawTower* tower = it->second;
0522           assert(tower);
0523
0524           const double energy_calib = tower->get_energy();
0525           hcalout_sum_energy_calib += energy_calib;
0526
0527           if(_is_sim) continue;
0528           const int col = tower->get_column();
0529           const int row = tower->get_row();
0530
0531           assert(
0532               _hcalout_recalib_const.find(make_pair(col, row)) != _hcalout_recalib_const.end());
0533           const double energy_recalib = energy_calib
0534               * _hcalout_recalib_const[make_pair(col, row)];
0535
0536               hcalout_sum_energy_recalib += energy_recalib;
0537          }
0538      }
0539
0540   const double EoP = emcal_sum_energy_recalib / abs(_eval_run.beam_mom);
0541   _eval_run.EoP = EoP;
0542
0543   // E/p
0544   if (_eval_run.good_e)
0545     {
0546       if (verbosity >= 3)
0547         cout << __PRETTY_FUNCTION__ << " sum_energy_calib = "
0548             << emcal_sum_energy_calib << " sum_energy_recalib = " << emcal_sum_energy_recalib
0549             << " _eval_run.beam_mom = " << _eval_run.beam_mom << endl;
0550
0551       TH2F * hEoP = dynamic_cast<TH2F *>(hm->getHisto("hEoP"));
0552       assert(hEoP);
0553
0554       hEoP->Fill(EoP, abs(_eval_run.beam_mom));
0555     }
0556
0557   // calibration file
0558    assert(fdata.is_open());
0559    fdata << sdata.str();
0560    fdata << endl;
0561
0562   _shower.emcal_e = emcal_sum_energy_calib;
0563   _shower.hcalin_e = hcalin_sum_energy_calib;
0564   _shower.hcalout_e = hcalout_sum_energy_calib;
0565   _shower.sum_e = emcal_sum_energy_calib + hcalin_sum_energy_calib + hcalout_sum_energy_calib ;
0566   _shower.hcal_asym = (hcalin_sum_energy_calib - hcalout_sum_energy_calib)/(hcalin_sum_energy_calib + hcalout_sum_energy_calib);
0567
0568   _shower.emcal_e_recal = emcal_sum_energy_recalib;
0569   _shower.hcalin_e_recal = hcalin_sum_energy_recalib;
0570   _shower.hcalout_e_recal = hcalout_sum_energy_recalib;
0571   _shower.sum_e_recal = emcal_sum_energy_recalib + hcalin_sum_energy_recalib + hcalout_sum_energy_recalib ;
0572
0573
0574   if(_fill_time_samples && !_is_sim)
0575   {
0576    //HCALIN
0577    {
0578     auto range = TOWER_RAW_LG_HCALIN->getTowers();
0579     for (auto it = range.first; it != range.second; ++it)
0580     {
0581       RawTower_Prototype2* tower = dynamic_cast<RawTower_Prototype2 *>(it->second);
0582       assert(tower);
0583
0584       int col = tower->get_column();
0585       int row = tower->get_row();
0586       int towerid = row + 4*col;
0587       for(int isample=0; isample<24; isample++)
0588        _time_samples.hcalin_time_samples[towerid][isample] =
0589                 tower->get_signal_samples(isample);
0590
0591      }
0592     }
0593
0594    //HCALOUT
0595    {
0596     auto range = TOWER_RAW_LG_HCALOUT->getTowers();
0597     for (auto it = range.first; it != range.second; ++it)
0598     {
0599       RawTower_Prototype2* tower = dynamic_cast<RawTower_Prototype2 *>(it->second);
0600       assert(tower);
0601
0602       int col = tower->get_column();
0603       int row = tower->get_row();
0604       int towerid = row + 4*col;
0605       for(int isample=0; isample<24; isample++)
0606        _time_samples.hcalout_time_samples[towerid][isample] =
0607         tower->get_signal_samples(isample);
0608
0609      }
0610     }
0611
0612    //EMCAL
0613    {
0614     auto range = TOWER_RAW_CEMC->getTowers();
0615     for (auto it = range.first; it != range.second; ++it)
0616     {
0617       RawTower_Prototype2* tower = dynamic_cast<RawTower_Prototype2 *>(it->second);
0618       assert(tower);
0619
0620       int col = tower->get_column();
0621       int row = tower->get_row();
0622       int towerid = row + 8*col;
0623       for(int isample=0; isample<24; isample++)
0624        _time_samples.emcal_time_samples[towerid][isample] =
0625                 tower->get_signal_samples(isample);
0626
0627      }
0628     }
0629   }
0630
0631   TTree * T = dynamic_cast<TTree *>(hm->getHisto("T"));
0632   assert(T);
0633   T->Fill();
0634
0635   return Fun4AllReturnCodes::EVENT_OK;
0636 }
0637
0638 int
0639 Proto2ShowerCalib::LoadRecalibMap(const std::string & file)
0640 {
0641   if (verbosity)
0642     {
0643       cout << __PRETTY_FUNCTION__ << " - input recalibration constant from "
0644           << file << endl;
0645     }
0646
0647   ifstream fcalib(file);
0648
0649   assert(fcalib.is_open());
0650
0651   string line;
0652
0653   while (not fcalib.eof())
0654     {
0655       getline(fcalib, line);
0656
0657       if (verbosity > 10)
0658         {
0659           cout << __PRETTY_FUNCTION__ << " get line " << line << endl;
0660         }
0661       istringstream sline(line);
0662
0663       int col = -1;
0664       int row = -1;
0665       double calib = 0;
0666
0667       sline >> col >> row >> calib;
0668
0669       if (not sline.fail())
0670         {
0671           if (verbosity)
0672             {
0673               cout << __PRETTY_FUNCTION__ << " - recalibration constant  "
0674                   << col << "," << row << " = " << calib << endl;
0675             }
0676
0677           _emcal_recalib_const[make_pair(col, row)] = calib;
0678         }
0679
0680     }
0681
0682   return _emcal_recalib_const.size();
0683 }
0684