HCAL/ShowerCalib/Proto3ShowerCalib.C

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