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