File indexing completed on 2025-08-05 08:11:13
0001 #include "EmulatorHistos.h"
0002
0003
0004 #include <calobase/TowerInfo.h>
0005 #include <calobase/TowerInfoContainer.h>
0006 #include <calobase/TowerInfoDefs.h>
0007
0008 #include <calotrigger/LL1Out.h>
0009 #include <calotrigger/LL1Outv1.h>
0010 #include <calotrigger/TriggerPrimitive.h>
0011 #include <calotrigger/TriggerPrimitivev1.h>
0012 #include <calotrigger/TriggerPrimitiveContainer.h>
0013 #include <calotrigger/TriggerPrimitiveContainerv1.h>
0014 #include <calotrigger/TriggerDefs.h>
0015
0016 #include <qautils/QAHistManagerDef.h>
0017
0018 #include <fun4all/Fun4AllHistoManager.h>
0019 #include <fun4all/Fun4AllReturnCodes.h>
0020
0021 #include <phool/getClass.h>
0022 #include <phool/phool.h> // for PHWHERE
0023
0024 #include <TFile.h>
0025 #include <TH1.h>
0026 #include <TH2.h>
0027 #include <TLorentzVector.h>
0028 #include <TNtuple.h>
0029 #include <TProfile2D.h>
0030 #include <TSystem.h>
0031 #include <TTree.h>
0032
0033 #include <cmath> // for log10, pow, sqrt, abs, M_PI
0034 #include <iostream> // for operator<<, endl, basic_...
0035 #include <limits>
0036 #include <map> // for operator!=, _Rb_tree_con...
0037 #include <string>
0038 #include <utility> // for pair
0039
0040 EmulatorHistos::EmulatorHistos(const std::string& name)
0041 : SubsysReco(name)
0042 {
0043 }
0044
0045 int EmulatorHistos::Init(PHCompositeNode* )
0046 {
0047 auto hm = QAHistManagerDef::getHistoManager();
0048 assert(hm);
0049
0050
0051 if (m_debug)
0052 {
0053 std::cout << "In EmulatorHistos::Init" << std::endl;
0054 }
0055
0056 TH1D *h = nullptr;
0057 TH2D *h2 = nullptr;
0058
0059
0060
0061 h = new TH1D("h_rejection_photon",";Emulated Threshold;Rejection Factor", 50, 0, 50);
0062 hm->registerHisto(h);
0063 h = new TH1D("h_rejection_jet",";Emulated Threshold;Rejection Factor", 50, 0, 50);
0064 hm->registerHisto(h);
0065
0066 h2 = new TH2D("h_photon_lutsum",";Max Digital Sum Photon;Corresponding 0.2x0.2 sum E [GeV]", 256, -0.5, 255.5, 50, 0, 50);
0067 hm->registerHisto(h2);
0068
0069 h2 = new TH2D("h_jet_lutsum",";Max Digital Sum Jet;Corresponding 0.8x0.8 sum E [GeV]", 256, -0.5, 255.5, 50, 0, 50);
0070 hm->registerHisto(h2);
0071
0072 for (int i = 0; i < 50; i++)
0073 {
0074 h = new TH1D(Form("h_eT_photon_%d", i),";Energy;", 50, 0, 50);
0075 hm->registerHisto(h);
0076 h = new TH1D(Form("h_eT_jet_%d", i),";Energy;", 50, 0, 50);
0077 hm->registerHisto(h);
0078 }
0079
0080 i_event = 0;
0081 h = new TH1D("h_nevent", "",1, 0, 2);
0082 hm->registerHisto(h);
0083
0084
0085 if (m_debug)
0086 {
0087 std::cout << "Leaving EmulatorHistos::Init" << std::endl;
0088 }
0089 return 0;
0090 }
0091
0092 int EmulatorHistos::process_event(PHCompositeNode* topNode)
0093 {
0094 i_event++;
0095
0096 return process_towers(topNode);
0097
0098 }
0099
0100 int EmulatorHistos::process_towers(PHCompositeNode* topNode)
0101 {
0102
0103
0104 TowerInfoContainer* towers_emcal = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_CEMC");
0105
0106 TowerInfoContainer* towers_hcalin = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_HCALIN");
0107
0108 TowerInfoContainer* towers_hcalout = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_HCALOUT");
0109
0110 LL1Out* ll1out_jet = findNode::getClass<LL1Out>(topNode, "LL1OUT_JET");
0111
0112 LL1Out* ll1out_photon = findNode::getClass<LL1Out>(topNode, "LL1OUT_PHOTON");
0113
0114
0115 auto hm = QAHistManagerDef::getHistoManager();
0116 assert(hm);
0117
0118 int ll1_photon[12][32];
0119 int ll1_jet[12][32];
0120 for (int i = 0; i < 12; i++)
0121 {
0122 for (int j = 0; j < 32; j++)
0123 {
0124 ll1_photon[i][j] = 0;
0125 ll1_jet[i][j] = 0;
0126 }
0127 }
0128 int max_photon = 0;
0129 int max_jet = 0;
0130
0131 if (ll1out_photon)
0132 {
0133 auto triggered_sums = ll1out_photon->getTriggeredSums();
0134
0135 for (auto &keypair : triggered_sums)
0136 {
0137 TriggerDefs::TriggerSumKey key = keypair.first;
0138 unsigned short bit = keypair.second;
0139 int phi = TriggerDefs::getSumPhiId(key) + TriggerDefs::getPrimitivePhiId_from_TriggerSumKey(key)*2;
0140 int eta = TriggerDefs::getSumEtaId(key);
0141 ll1_photon[eta][phi] = static_cast<int>(bit);
0142 if (static_cast<int>(bit) > max_photon)
0143 max_photon = static_cast<int>(bit);
0144 }
0145 }
0146 if (ll1out_jet)
0147 {
0148 auto triggered_sums = ll1out_jet->getTriggeredSums();
0149 for (auto &keypair : triggered_sums)
0150 {
0151 TriggerDefs::TriggerSumKey key = keypair.first;
0152 unsigned short bit = keypair.second;
0153 int phi = static_cast<int>(key & 0xffffU);
0154 int eta = static_cast<int>((key >> 16U) & 0xffffU);
0155 ll1_jet[eta][phi] = static_cast<int>(bit);
0156 if (static_cast<int>(bit) > max_jet)
0157 max_jet = static_cast<int>(bit);
0158 }
0159 }
0160
0161 float emcal_energies[12][35]{0.0};
0162 float hcal_energies[12][35]{0.0};
0163 float max_energy_emcal = 0.0;
0164 int size;
0165 if (towers_emcal)
0166 {
0167 size = towers_emcal->size();
0168 for (int itower = 0; itower < size ; itower++)
0169 {
0170 TowerInfo* tower = towers_emcal->get_tower_at_channel(itower);
0171
0172 float offlineenergy = tower->get_energy();
0173 unsigned int towerkey = towers_emcal->encode_key(itower);
0174 int ieta = towers_emcal->getTowerEtaBin(towerkey);
0175 int iphi = towers_emcal->getTowerPhiBin(towerkey);
0176 emcal_energies[ieta/8][iphi/8] += offlineenergy;
0177 }
0178 for (int i = 0; i < 12; i++)
0179 {
0180 for (int j = 0; j < 32; j++)
0181 {
0182 if (emcal_energies[i][j] > max_energy_emcal)
0183 {
0184 max_energy_emcal = emcal_energies[i][j];
0185 }
0186 }
0187 }
0188 }
0189
0190 if (towers_hcalin)
0191 {
0192 size = towers_hcalin->size();
0193 for (int itower = 0; itower < size ; itower++)
0194 {
0195 TowerInfo* tower = towers_hcalin->get_tower_at_channel(itower);
0196
0197 float offlineenergy = tower->get_energy();
0198 unsigned int towerkey = towers_hcalin->encode_key(itower);
0199 int ieta = towers_hcalin->getTowerEtaBin(towerkey);
0200 int iphi = towers_hcalin->getTowerPhiBin(towerkey);
0201 hcal_energies[ieta/2][iphi/2] += offlineenergy;
0202 }
0203
0204 }
0205 if (towers_hcalout)
0206 {
0207 size = towers_hcalout->size();
0208 for (int itower = 0; itower < size ; itower++)
0209 {
0210 TowerInfo* tower = towers_hcalout->get_tower_at_channel(itower);
0211
0212 float offlineenergy = tower->get_energy();
0213 unsigned int towerkey = towers_hcalout->encode_key(itower);
0214 int ieta = towers_hcalout->getTowerEtaBin(towerkey);
0215 int iphi = towers_hcalout->getTowerPhiBin(towerkey);
0216 hcal_energies[ieta/2][iphi/2] += offlineenergy;
0217 }
0218
0219 }
0220
0221
0222 float jet_energies[9][32]{};
0223 float max_energy_jetpatch = 0.0;
0224
0225 for (int i = 0; i < 3; i++)
0226 {
0227 for (int j = 0; j < 12; j++)
0228 {
0229 emcal_energies[j][i + 32] = emcal_energies[j][i];
0230 hcal_energies[j][i + 32] = hcal_energies[j][i];
0231 }
0232 }
0233 for (int i = 0; i < 9; i++)
0234 {
0235 for (int j = 0; j < 32; j++)
0236 {
0237 for (int k = 0; k < 16; k++)
0238 {
0239 jet_energies[i][j] += emcal_energies[i + k % 4][j + k / 4];
0240 jet_energies[i][j] += hcal_energies[i + k % 4][j + k / 4];
0241 }
0242 if (jet_energies[i][j] > max_energy_jetpatch)
0243 {
0244 max_energy_jetpatch = jet_energies[i][j];
0245 }
0246 }
0247 }
0248
0249 auto h_rej_photon = dynamic_cast<TH1*>(hm->getHisto("h_rejection_photon"));
0250 for (int i = 0; i < max_photon; i++)
0251 {
0252 h_rej_photon->Fill(i);
0253 }
0254
0255 auto h_rej_jet = dynamic_cast<TH1*>(hm->getHisto("h_rejection_jet"));
0256 for (int i = 0; i < max_jet; i++)
0257 {
0258 h_rej_jet->Fill(i);
0259 }
0260
0261 auto h_photon_lutsum = dynamic_cast<TH2*>(hm->getHisto("h_photon_lutsum"));
0262 auto h_jet_lutsum = dynamic_cast<TH2*>(hm->getHisto("h_jet_lutsum"));
0263 for (int i = 0; i < 12; i++)
0264 {
0265 for (int j = 0; j < 32; j++)
0266 {
0267 h_photon_lutsum->Fill(ll1_photon[i][j], emcal_energies[i][j]);
0268 if (i < 9)
0269 {
0270 h_jet_lutsum->Fill(ll1_jet[i][j], jet_energies[i][j]);
0271 }
0272 }
0273 }
0274
0275
0276 for ( int i = 0; i < 50; i++)
0277 {
0278 auto h_eT_photon = dynamic_cast<TH1*>(hm->getHisto(Form("h_eT_photon_%d", i)));
0279 if (max_photon >= i)
0280 h_eT_photon->Fill(max_energy_emcal);
0281 auto h_eT_jet = dynamic_cast<TH1*>(hm->getHisto(Form("h_eT_jet_%d", i)));
0282 if (max_jet >= i)
0283 h_eT_jet->Fill(max_energy_emcal);
0284 }
0285 auto h_nevent = dynamic_cast<TH1*>(hm->getHisto("h_nevent"));
0286 h_nevent->Fill(1);
0287 return Fun4AllReturnCodes::EVENT_OK;
0288 }
0289
0290 int EmulatorHistos::End(PHCompositeNode*)
0291 {
0292
0293
0294
0295
0296
0297
0298
0299
0300
0301
0302
0303
0304
0305
0306
0307
0308
0309
0310
0311
0312
0313
0314
0315
0316
0317
0318
0319 return Fun4AllReturnCodes::EVENT_OK;
0320 }
