File indexing completed on 2025-08-05 08:14:36
0001 #include "Proto2ShowerCalib.h"
0002 #include "TemperatureCorrection.h"
0003
0004 #include <prototype2/RawTower_Temperature.h>
0005 #include <prototype2/RawTower_Prototype2.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(Proto2ShowerCalib::Eval_Cluster);
0045 ClassImp(Proto2ShowerCalib::Eval_Run);
0046
0047 Proto2ShowerCalib::Proto2ShowerCalib(const std::string &filename) :
0048 SubsysReco("Proto2ShowerCalib"), _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 Proto2ShowerCalib::~Proto2ShowerCalib()
0073 {
0074 }
0075
0076 Fun4AllHistoManager *
0077 Proto2ShowerCalib::get_HistoManager()
0078 {
0079
0080 Fun4AllServer *se = Fun4AllServer::instance();
0081 Fun4AllHistoManager *hm = se->getHistoManager("Proto2ShowerCalib_HISTOS");
0082
0083 if (not hm)
0084 {
0085 cout
0086 << "Proto2ShowerCalib::get_HistoManager - Making Fun4AllHistoManager Proto2ShowerCalib_HISTOS"
0087 << endl;
0088 hm = new Fun4AllHistoManager("Proto2ShowerCalib_HISTOS");
0089 se->registerHistoManager(hm);
0090 }
0091
0092 assert(hm);
0093
0094 return hm;
0095 }
0096
0097 int
0098 Proto2ShowerCalib::InitRun(PHCompositeNode *topNode)
0099 {
0100 if (verbosity)
0101 cout << "Proto2ShowerCalib::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 Proto2ShowerCalib::End(PHCompositeNode *topNode)
0120 {
0121 cout << "Proto2ShowerCalib::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
0130
0131
0132 fdata.close();
0133
0134 return Fun4AllReturnCodes::EVENT_OK;
0135 }
0136
0137 int
0138 Proto2ShowerCalib::Init(PHCompositeNode *topNode)
0139 {
0140
0141 _ievent = 0;
0142
0143 cout << "Proto2ShowerCalib::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
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 Proto2ShowerCalib::process_event(PHCompositeNode *topNode)
0203 {
0204
0205 if (verbosity > 2)
0206 cout << "Proto2ShowerCalib::process_event() entered" << endl;
0207
0208
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
0240 EventHeader* eventheader = findNode::getClass<EventHeader>(topNode,
0241 "EventHeader");
0242 if (not _is_sim)
0243 {
0244 assert(eventheader);
0245
0246 _eval_run.run = eventheader->get_RunNumber();
0247 if (verbosity > 4)
0248 cout << __PRETTY_FUNCTION__ << _eval_run.run << endl;
0249
0250 _eval_run.event = eventheader->get_EvtSequence();
0251 }
0252
0253 if (_is_sim)
0254 {
0255
0256 PHG4TruthInfoContainer* truthInfoList = findNode::getClass<
0257 PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
0258
0259 assert(truthInfoList);
0260
0261 _eval_run.run = -1;
0262
0263 const PHG4Particle * p = truthInfoList->GetPrimaryParticleRange().first->second;
0264 assert(p);
0265
0266 const PHG4VtxPoint * v = truthInfoList->GetVtx(p->get_vtx_id());
0267 assert(v);
0268
0269 _eval_run.beam_mom = sqrt(
0270 p->get_px() * p->get_px() + p->get_py() * p->get_py()
0271 + p->get_pz() * p->get_pz());
0272 _eval_run.truth_y = v->get_y();
0273 _eval_run.truth_z = v->get_z();
0274
0275 }
0276
0277
0278 TH1F * hNormalization = dynamic_cast<TH1F *>(hm->getHisto("hNormalization"));
0279 assert(hNormalization);
0280
0281 hNormalization->Fill("ALL", 1);
0282
0283 RawTowerContainer* TOWER_RAW_CEMC = findNode::getClass<RawTowerContainer>(
0284 topNode, _is_sim ? "TOWER_RAW_LG_CEMC" : "TOWER_RAW_CEMC");
0285 assert(TOWER_RAW_CEMC);
0286 RawTowerContainer* TOWER_CALIB_CEMC = findNode::getClass<RawTowerContainer>(
0287 topNode, _is_sim ? "TOWER_CALIB_LG_CEMC" : "TOWER_CALIB_CEMC");
0288 assert(TOWER_CALIB_CEMC);
0289
0290
0291 RawTowerContainer* TOWER_CALIB_TRIGGER_VETO = findNode::getClass<
0292 RawTowerContainer>(topNode, "TOWER_CALIB_TRIGGER_VETO");
0293 if (not _is_sim)
0294 {
0295 assert(TOWER_CALIB_TRIGGER_VETO);
0296 }
0297
0298 RawTowerContainer* TOWER_CALIB_HODO_HORIZONTAL = findNode::getClass<
0299 RawTowerContainer>(topNode, "TOWER_CALIB_HODO_HORIZONTAL");
0300 if (not _is_sim)
0301 {
0302 assert(TOWER_CALIB_HODO_HORIZONTAL);
0303 }
0304 RawTowerContainer* TOWER_CALIB_HODO_VERTICAL = findNode::getClass<
0305 RawTowerContainer>(topNode, "TOWER_CALIB_HODO_VERTICAL");
0306 if (not _is_sim)
0307 {
0308 assert(TOWER_CALIB_HODO_VERTICAL);
0309 }
0310
0311 RawTowerContainer* TOWER_TEMPERATURE_EMCAL = findNode::getClass<
0312 RawTowerContainer>(topNode, "TOWER_TEMPERATURE_EMCAL");
0313 if (not _is_sim)
0314 {
0315 assert(TOWER_TEMPERATURE_EMCAL);
0316 }
0317
0318 RawTowerContainer* TOWER_CALIB_C1 = findNode::getClass<RawTowerContainer>(
0319 topNode, "TOWER_CALIB_C1");
0320 if (not _is_sim)
0321 {
0322 assert(TOWER_CALIB_C1);
0323 }
0324 RawTowerContainer* TOWER_CALIB_C2 = findNode::getClass<RawTowerContainer>(
0325 topNode, "TOWER_CALIB_C2");
0326 if (not _is_sim)
0327 {
0328 assert(TOWER_CALIB_C2);
0329 }
0330
0331
0332 bool cherekov_e = false;
0333 if (not _is_sim)
0334 {
0335 RawTower * t_c2_in = NULL;
0336 RawTower * t_c2_out = NULL;
0337
0338 assert(eventheader);
0339 if (eventheader->get_RunNumber() >= 2105)
0340 {
0341 t_c2_in = TOWER_CALIB_C2->getTower(10);
0342 t_c2_out = TOWER_CALIB_C2->getTower(11);
0343 }
0344 else
0345 {
0346 t_c2_in = TOWER_CALIB_C2->getTower(0);
0347 t_c2_out = TOWER_CALIB_C2->getTower(1);
0348 }
0349 assert(t_c2_in);
0350 assert(t_c2_out);
0351
0352 const double c2_in = t_c2_in->get_energy();
0353 const double c2_out = t_c2_out->get_energy();
0354 const double c1 = TOWER_CALIB_C1->getTower(0)->get_energy();
0355
0356 _eval_run.C2_sum = c2_in + c2_out;
0357 cherekov_e = (_eval_run.C2_sum) > 100;
0358 hNormalization->Fill("C2-e", cherekov_e);
0359
0360 TH2F * hCheck_Cherenkov = dynamic_cast<TH2F *>(hm->getHisto(
0361 "hCheck_Cherenkov"));
0362 assert(hCheck_Cherenkov);
0363 hCheck_Cherenkov->Fill(c1, "C1", 1);
0364 hCheck_Cherenkov->Fill(c2_in, "C2 in", 1);
0365 hCheck_Cherenkov->Fill(c2_out, "C2 out", 1);
0366 hCheck_Cherenkov->Fill(c2_in + c2_out, "C2 sum", 1);
0367 }
0368
0369
0370 TH1F * hCheck_Veto = dynamic_cast<TH1F *>(hm->getHisto("hCheck_Veto"));
0371 assert(hCheck_Veto);
0372 bool trigger_veto_pass = true;
0373 if (not _is_sim)
0374 {
0375 auto range = TOWER_CALIB_TRIGGER_VETO->getTowers();
0376 for (auto it = range.first; it != range.second; ++it)
0377 {
0378 RawTower* tower = it->second;
0379 assert(tower);
0380
0381 hCheck_Veto->Fill(tower->get_energy());
0382
0383 if (abs(tower->get_energy()) > 15)
0384 trigger_veto_pass = false;
0385 }
0386 }
0387 hNormalization->Fill("trigger_veto_pass", trigger_veto_pass);
0388 _eval_run.trigger_veto_pass = trigger_veto_pass;
0389
0390
0391 TH1F * hCheck_Hodo_H = dynamic_cast<TH1F *>(hm->getHisto("hCheck_Hodo_H"));
0392 assert(hCheck_Hodo_H);
0393 int hodo_h_count = 0;
0394 if (not _is_sim)
0395 {
0396 auto range = TOWER_CALIB_HODO_HORIZONTAL->getTowers();
0397 for (auto it = range.first; it != range.second; ++it)
0398 {
0399 RawTower* tower = it->second;
0400 assert(tower);
0401
0402 hCheck_Hodo_H->Fill(tower->get_energy());
0403
0404 if (abs(tower->get_energy()) > 30)
0405 {
0406 hodo_h_count++;
0407 _eval_run.hodo_h = tower->get_id();
0408 }
0409 }
0410 }
0411 const bool valid_hodo_h = hodo_h_count == 1;
0412 hNormalization->Fill("valid_hodo_h", valid_hodo_h);
0413 _eval_run.valid_hodo_h = valid_hodo_h;
0414
0415 TH1F * hCheck_Hodo_V = dynamic_cast<TH1F *>(hm->getHisto("hCheck_Hodo_V"));
0416 assert(hCheck_Hodo_V);
0417 int hodo_v_count = 0;
0418 if (not _is_sim)
0419 {
0420 auto range = TOWER_CALIB_HODO_VERTICAL->getTowers();
0421 for (auto it = range.first; it != range.second; ++it)
0422 {
0423 RawTower* tower = it->second;
0424 assert(tower);
0425
0426 hCheck_Hodo_V->Fill(tower->get_energy());
0427
0428 if (abs(tower->get_energy()) > 30)
0429 {
0430 hodo_v_count++;
0431 _eval_run.hodo_v = tower->get_id();
0432 }
0433 }
0434 }
0435 const bool valid_hodo_v = hodo_v_count == 1;
0436 _eval_run.valid_hodo_v = valid_hodo_v;
0437 hNormalization->Fill("valid_hodo_v", valid_hodo_v);
0438
0439 const bool good_e = (valid_hodo_v and valid_hodo_h and cherekov_e
0440 and trigger_veto_pass) and (not _is_sim);
0441 hNormalization->Fill("good_e", good_e);
0442
0443
0444 pair<int, int> max_3x3 = find_max(TOWER_CALIB_CEMC, 3);
0445 pair<int, int> max_5x5 = find_max(TOWER_CALIB_CEMC, 5);
0446
0447 _eval_3x3_raw.max_col = max_3x3.first;
0448 _eval_3x3_raw.max_row = max_3x3.second;
0449 _eval_3x3_prod.max_col = max_3x3.first;
0450 _eval_3x3_prod.max_row = max_3x3.second;
0451 _eval_3x3_temp.max_col = max_3x3.first;
0452 _eval_3x3_temp.max_row = max_3x3.second;
0453 _eval_3x3_recalib.max_col = max_3x3.first;
0454 _eval_3x3_recalib.max_row = max_3x3.second;
0455
0456 _eval_5x5_raw.max_col = max_5x5.first;
0457 _eval_5x5_raw.max_row = max_5x5.second;
0458 _eval_5x5_prod.max_col = max_5x5.first;
0459 _eval_5x5_prod.max_row = max_5x5.second;
0460 _eval_5x5_temp.max_col = max_5x5.first;
0461 _eval_5x5_temp.max_row = max_5x5.second;
0462 _eval_5x5_recalib.max_col = max_5x5.first;
0463 _eval_5x5_recalib.max_row = max_5x5.second;
0464
0465
0466 bool good_temp = true;
0467 double sum_energy_calib = 0;
0468 double sum_energy_T = 0;
0469 TH1F * hTemperature = dynamic_cast<TH1F *>(hm->getHisto("hTemperature"));
0470 assert(hTemperature);
0471
0472 stringstream sdata;
0473
0474 if (good_e)
0475 sdata << abs(_eval_run.beam_mom) << "\t";
0476
0477
0478 {
0479 auto range = TOWER_CALIB_CEMC->getTowers();
0480 for (auto it = range.first; it != range.second; ++it)
0481 {
0482
0483 RawTowerDefs::keytype key = it->first;
0484 RawTower* tower = it->second;
0485 assert(tower);
0486
0487 const int col = tower->get_bineta();
0488 const int row = tower->get_binphi();
0489
0490 if (col < 0 or col >= 8)
0491 continue;
0492 if (row < 0 or row >= 8)
0493 continue;
0494
0495 const double energy_calib = tower->get_energy();
0496 sum_energy_calib += energy_calib;
0497
0498 RawTower* tower_raw = TOWER_RAW_CEMC->getTower(key);
0499 assert(tower_raw);
0500
0501 double energy_T = 0;
0502 if (not _is_sim)
0503 {
0504 RawTower_Temperature * temp_t =
0505 dynamic_cast<RawTower_Temperature *>(TOWER_TEMPERATURE_EMCAL->getTower(
0506 tower->get_row(), tower->get_column()));
0507 assert(temp_t);
0508
0509 const double T = temp_t->get_temperature_from_time(
0510 eventheader->get_TimeStamp());
0511 hTemperature->Fill(T);
0512
0513 if (T < 25 or T > 35)
0514 good_temp = false;
0515
0516 energy_T = TemperatureCorrection::Apply(energy_calib, T);
0517 }
0518
0519
0520 assert(
0521 _recalib_const.find(make_pair(col, row)) != _recalib_const.end());
0522 const double energy_recalib = energy_T
0523 * _recalib_const[make_pair(col, row)];
0524
0525
0526 sum_energy_T += energy_T;
0527
0528
0529
0530
0531 if (col >= max_5x5.first - 2 and col <= max_5x5.first + 2
0532 and row >= max_5x5.second - 2 and row <= max_5x5.second + 2)
0533 {
0534 sdata << tower->get_energy() << "\t";
0535 }
0536 else
0537 {
0538 sdata << 0 << "\t";
0539 }
0540
0541
0542 if (col >= max_3x3.first - 1 and col <= max_3x3.first + 1)
0543 if (row >= max_3x3.second - 1 and row <= max_3x3.second + 1)
0544 {
0545
0546
0547 _eval_3x3_raw.average_col += abs(tower_raw->get_energy()) * col;
0548 _eval_3x3_raw.average_row += abs(tower_raw->get_energy()) * row;
0549 _eval_3x3_raw.sum_E += abs(tower_raw->get_energy());
0550
0551 _eval_3x3_prod.average_col += energy_calib * col;
0552 _eval_3x3_prod.average_row += energy_calib * row;
0553 _eval_3x3_prod.sum_E += energy_calib;
0554
0555 _eval_3x3_temp.average_col += energy_T * col;
0556 _eval_3x3_temp.average_row += energy_T * row;
0557 _eval_3x3_temp.sum_E += energy_T;
0558
0559 _eval_3x3_recalib.average_col += energy_recalib * col;
0560 _eval_3x3_recalib.average_row += energy_recalib * row;
0561 _eval_3x3_recalib.sum_E += energy_recalib;
0562 }
0563
0564
0565 if (col >= max_5x5.first - 2 and col <= max_5x5.first + 2)
0566 if (row >= max_5x5.second - 2 and row <= max_5x5.second + 2)
0567 {
0568
0569
0570 _eval_5x5_raw.average_col += abs(tower_raw->get_energy()) * col;
0571 _eval_5x5_raw.average_row += abs(tower_raw->get_energy()) * row;
0572 _eval_5x5_raw.sum_E += abs(tower_raw->get_energy());
0573
0574 _eval_5x5_prod.average_col += energy_calib * col;
0575 _eval_5x5_prod.average_row += energy_calib * row;
0576 _eval_5x5_prod.sum_E += energy_calib;
0577
0578 _eval_5x5_temp.average_col += energy_T * col;
0579 _eval_5x5_temp.average_row += energy_T * row;
0580 _eval_5x5_temp.sum_E += energy_T;
0581
0582 _eval_5x5_recalib.average_col += energy_recalib * col;
0583 _eval_5x5_recalib.average_row += energy_recalib * row;
0584 _eval_5x5_recalib.sum_E += energy_recalib;
0585 }
0586 }
0587 }
0588
0589 _eval_3x3_raw.reweight_clus_pol();
0590 _eval_5x5_raw.reweight_clus_pol();
0591 _eval_3x3_prod.reweight_clus_pol();
0592 _eval_5x5_prod.reweight_clus_pol();
0593 _eval_3x3_temp.reweight_clus_pol();
0594 _eval_5x5_temp.reweight_clus_pol();
0595 _eval_3x3_recalib.reweight_clus_pol();
0596 _eval_5x5_recalib.reweight_clus_pol();
0597
0598 const double EoP = sum_energy_T / abs(_eval_run.beam_mom);
0599 hNormalization->Fill("good_temp", good_temp);
0600
0601 bool good_data = good_e and good_temp;
0602 hNormalization->Fill("good_data", good_data);
0603
0604 _eval_run.good_temp = good_temp;
0605 _eval_run.good_e = good_e;
0606 _eval_run.good_data = good_data;
0607 _eval_run.sum_energy_T = sum_energy_T;
0608 _eval_run.EoP = EoP;
0609
0610
0611 if (good_data)
0612 {
0613 if (verbosity >= 3)
0614 cout << __PRETTY_FUNCTION__ << " sum_energy_calib = "
0615 << sum_energy_calib << " sum_energy_T = " << sum_energy_T
0616 << " _eval_run.beam_mom = " << _eval_run.beam_mom << endl;
0617
0618 TH2F * hEoP = dynamic_cast<TH2F *>(hm->getHisto("hEoP"));
0619 assert(hEoP);
0620
0621 hEoP->Fill(EoP, abs(_eval_run.beam_mom));
0622 }
0623
0624
0625 if (good_data and abs(_eval_run.beam_mom) >= 4
0626 and abs(_eval_run.beam_mom) <= 8)
0627 {
0628 assert(fdata.is_open());
0629
0630 fdata << sdata.str();
0631
0632 fdata << endl;
0633 }
0634
0635 TTree * T = dynamic_cast<TTree *>(hm->getHisto("T"));
0636 assert(T);
0637 T->Fill();
0638
0639 return Fun4AllReturnCodes::EVENT_OK;
0640 }
0641
0642 pair<int, int>
0643 Proto2ShowerCalib::find_max(RawTowerContainer* towers, int cluster_size)
0644 {
0645 const int clus_edge_size = (cluster_size - 1) / 2;
0646 assert(clus_edge_size >= 0);
0647
0648 pair<int, int> max(-1000, -1000);
0649 double max_e = 0;
0650
0651 for (int col = 0; col < n_size; ++col)
0652 for (int row = 0; row < n_size; ++row)
0653 {
0654 double energy = 0;
0655
0656 for (int dcol = col - clus_edge_size; dcol <= col + clus_edge_size;
0657 ++dcol)
0658 for (int drow = row - clus_edge_size; drow <= row + clus_edge_size;
0659 ++drow)
0660 {
0661 if (dcol < 0 or drow < 0)
0662 continue;
0663
0664 RawTower * t = towers->getTower(dcol, drow);
0665 if (t)
0666 energy += t->get_energy();
0667 }
0668
0669 if (energy > max_e)
0670 {
0671 max = make_pair(col, row);
0672 max_e = energy;
0673 }
0674 }
0675
0676 return max;
0677 }
0678
0679 int
0680 Proto2ShowerCalib::LoadRecalibMap(const std::string & file)
0681 {
0682 if (verbosity)
0683 {
0684 cout << __PRETTY_FUNCTION__ << " - input recalibration constant from "
0685 << file << endl;
0686 }
0687
0688 ifstream fcalib(file);
0689
0690 assert(fcalib.is_open());
0691
0692 string line;
0693
0694 while (not fcalib.eof())
0695 {
0696 getline(fcalib, line);
0697
0698 if (verbosity > 10)
0699 {
0700 cout << __PRETTY_FUNCTION__ << " get line " << line << endl;
0701 }
0702 istringstream sline(line);
0703
0704 int col = -1;
0705 int row = -1;
0706 double calib = 0;
0707
0708 sline >> col >> row >> calib;
0709
0710 if (not sline.fail())
0711 {
0712 if (verbosity)
0713 {
0714 cout << __PRETTY_FUNCTION__ << " - recalibration constant "
0715 << col << "," << row << " = " << calib << endl;
0716 }
0717
0718 _recalib_const[make_pair(col, row)] = calib;
0719 }
0720
0721 }
0722
0723 return _recalib_const.size();
0724 }
0725
