EMCal/ShowerCalib/Proto3ShowerCalib.C

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