File indexing completed on 2025-08-06 08:14:15
0001
0002
0003
0004
0005
0006 #include "HCalJetPhiShift.h"
0007
0008 #include "g4eval/CaloEvalStack.h"
0009 #include "g4eval/CaloRawClusterEval.h"
0010 #include "g4eval/CaloRawTowerEval.h"
0011 #include "g4eval/CaloTruthEval.h"
0012
0013 #include <g4main/PHG4Particle.h>
0014 #include <g4main/PHG4Shower.h>
0015 #include <g4main/PHG4TruthInfoContainer.h>
0016 #include <g4main/PHG4VtxPoint.h>
0017
0018
0019
0020
0021
0022
0023
0024 #include <fun4all/Fun4AllReturnCodes.h>
0025 #include <fun4all/PHTFileServer.h>
0026
0027 #include <phool/PHCompositeNode.h>
0028 #include <phool/getClass.h>
0029
0030 #include <calobase/RawTower.h>
0031 #include <calobase/RawTowerContainer.h>
0032 #include <calobase/RawTowerGeom.h>
0033 #include <calobase/RawTowerGeomContainer.h>
0034 #include <calobase/TowerInfoContainer.h>
0035 #include <calobase/TowerInfo.h>
0036
0037 #include <TFile.h>
0038 #include <TTree.h>
0039
0040
0041 HCalJetPhiShift::HCalJetPhiShift(const std::string &name, const std::string &outputFile):
0042 SubsysReco(name),
0043 m_outputFileName(outputFile),
0044 m_T(nullptr),
0045 m_event(-1),
0046 m_nTow_in(0),
0047 m_nTow_out(0),
0048 m_nTow_emc(0),
0049 m_eta(),
0050 m_phi(),
0051 m_e(),
0052 m_pt(),
0053 m_vx(),
0054 m_vy(),
0055 m_vz(),
0056 m_id(),
0057 m_eta_in(),
0058 m_phi_in(),
0059 m_e_in(),
0060 m_ieta_in(),
0061 m_iphi_in(),
0062 m_eta_out(),
0063 m_phi_out(),
0064 m_e_out(),
0065 m_ieta_out(),
0066 m_iphi_out(),
0067 m_eta_emc(),
0068 m_phi_emc(),
0069 m_e_emc(),
0070 m_ieta_emc(),
0071 m_iphi_emc()
0072
0073 {
0074 std::cout << "HCalJetPhiShift::HCalJetPhiShift(const std::string &name) Calling ctor" << std::endl;
0075 }
0076
0077
0078 HCalJetPhiShift::~HCalJetPhiShift()
0079 {
0080 std::cout << "HCalJetPhiShift::~HCalJetPhiShift() Calling dtor" << std::endl;
0081 }
0082
0083
0084 int HCalJetPhiShift::Init(PHCompositeNode* )
0085 {
0086 std::cout << "HCalJetPhiShift::Init(PHCompositeNode *topNode) Initializing" << std::endl;
0087 PHTFileServer::get().open(m_outputFileName, "RECREATE");
0088 std::cout << "HCalJetPhiShift::Init - Output to " << m_outputFileName << std::endl;
0089
0090
0091 m_T = new TTree("T", "HCalJetPhiShift Tree");
0092 m_T->Branch("event", &m_event);
0093 m_T->Branch("eta", &m_eta);
0094 m_T->Branch("phi", &m_phi);
0095 m_T->Branch("e", &m_e);
0096 m_T->Branch("pt", &m_pt);
0097 m_T->Branch("vx", &m_vx);
0098 m_T->Branch("vy", &m_vy);
0099 m_T->Branch("vz", &m_vz);
0100 m_T->Branch("id", &m_id);
0101 m_T->Branch("nTow_in", &m_nTow_in);
0102 m_T->Branch("eta_in", &m_eta_in);
0103 m_T->Branch("phi_in", &m_phi_in);
0104 m_T->Branch("e_in", &m_e_in);
0105 m_T->Branch("ieta_in", &m_ieta_in);
0106 m_T->Branch("iphi_in", &m_iphi_in);
0107 m_T->Branch("nTow_out", &m_nTow_out);
0108 m_T->Branch("eta_out", &m_eta_out);
0109 m_T->Branch("phi_out", &m_phi_out);
0110 m_T->Branch("e_out", &m_e_out);
0111 m_T->Branch("ieta_out", &m_ieta_out);
0112 m_T->Branch("iphi_out", &m_iphi_out);
0113 m_T->Branch("nTow_emc", &m_nTow_emc);
0114 m_T->Branch("eta_emc", &m_eta_emc);
0115 m_T->Branch("phi_emc", &m_phi_emc);
0116 m_T->Branch("e_emc", &m_e_emc);
0117 m_T->Branch("ieta_emc", &m_ieta_emc);
0118 m_T->Branch("iphi_emc", &m_iphi_emc);
0119
0120
0121 return Fun4AllReturnCodes::EVENT_OK;
0122 }
0123
0124
0125 int HCalJetPhiShift::InitRun(PHCompositeNode *topNode)
0126 {
0127 std::cout << "HCalJetPhiShift::InitRun(PHCompositeNode *topNode) Initializing for Run XXX" << std::endl;
0128 return Fun4AllReturnCodes::EVENT_OK;
0129 }
0130
0131
0132 int HCalJetPhiShift::process_event(PHCompositeNode *topNode)
0133 {
0134
0135 ++m_event;
0136
0137
0138 FillTTree(topNode);
0139
0140 return Fun4AllReturnCodes::EVENT_OK;
0141 }
0142
0143
0144 int HCalJetPhiShift::ResetEvent(PHCompositeNode *topNode)
0145 {
0146
0147 m_nTow_in = 0;
0148 m_nTow_out = 0;
0149 m_nTow_emc = 0;
0150 m_id.clear();
0151 m_eta_in.clear();
0152 m_phi_in.clear();
0153 m_e_in.clear();
0154 m_ieta_in.clear();
0155 m_iphi_in.clear();
0156
0157 m_eta_out.clear();
0158 m_phi_out.clear();
0159 m_e_out.clear();
0160 m_ieta_out.clear();
0161 m_iphi_out.clear();
0162
0163 m_eta_emc.clear();
0164 m_phi_emc.clear();
0165 m_e_emc.clear();
0166 m_ieta_emc.clear();
0167 m_iphi_emc.clear();
0168
0169 return Fun4AllReturnCodes::EVENT_OK;
0170 }
0171
0172
0173 int HCalJetPhiShift::EndRun(const int runnumber)
0174 {
0175 std::cout << "HCalJetPhiShift::EndRun(const int runnumber) Ending Run for Run " << runnumber << std::endl;
0176 return Fun4AllReturnCodes::EVENT_OK;
0177 }
0178
0179
0180 int HCalJetPhiShift::End(PHCompositeNode *topNode)
0181 {
0182 std::cout << "HCalJetPhiShift::End - Output to " << m_outputFileName << std::endl;
0183 PHTFileServer::get().cd(m_outputFileName);
0184
0185 m_T->Write();
0186 std::cout << "HCalJetPhiShift::End(PHCompositeNode *topNode) This is the End..." << std::endl;
0187 return Fun4AllReturnCodes::EVENT_OK;
0188 }
0189
0190
0191 int HCalJetPhiShift::Reset(PHCompositeNode *topNode)
0192 {
0193 std::cout << "HCalJetPhiShift::Reset(PHCompositeNode *topNode) being Reset" << std::endl;
0194 return Fun4AllReturnCodes::EVENT_OK;
0195 }
0196
0197
0198 void HCalJetPhiShift::Print(const std::string &what) const
0199 {
0200 std::cout << "HCalJetPhiShift::Print(const std::string &what) const Printing info for " << what << std::endl;
0201 }
0202
0203
0204 int HCalJetPhiShift::FillTTree(PHCompositeNode *topNode)
0205 {
0206
0207
0208
0209 PHG4TruthInfoContainer* truthinfo = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
0210 if (!truthinfo)
0211 {
0212 std::cout << PHWHERE << " ERROR: Can't find G4TruthInfo" << std::endl;
0213 exit(-1);
0214 }
0215
0216 PHG4TruthInfoContainer::ConstRange range = truthinfo->GetPrimaryParticleRange();
0217 for (PHG4TruthInfoContainer::ConstIterator iter = range.first;
0218 iter != range.second;
0219 ++iter)
0220 {
0221 PHG4Particle* primary = iter->second;
0222
0223 if (primary->get_e() < 0.)
0224 {
0225 continue;
0226 }
0227
0228 float gpx = primary->get_px();
0229 float gpy = primary->get_py();
0230 float gpz = primary->get_pz();
0231 m_e = primary->get_e();
0232
0233 m_pt = std::sqrt(gpx * gpx + gpy * gpy);
0234 m_eta = NAN;
0235 if (m_pt != 0.0)
0236 {
0237 m_eta = std::asinh(gpz / m_pt);
0238 }
0239 m_phi = std::atan2(gpy, gpx);
0240
0241 PHG4VtxPoint* gvertex = truthinfo->GetPrimaryVtx(truthinfo->GetPrimaryVertexIndex());
0242 m_vx = gvertex->get_x();
0243 m_vy = gvertex->get_y();
0244 m_vz = gvertex->get_z();
0245 }
0246
0247
0248 TowerInfoContainer *towersIH3 = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_HCALIN");
0249 TowerInfoContainer *towersOH3 = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_HCALOUT");
0250
0251 TowerInfoContainer *towersEM3 = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_RAW_CEMC");
0252 RawTowerGeomContainer *tower_geomIH = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALIN");
0253 RawTowerGeomContainer *tower_geomOH = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALOUT");
0254 RawTowerGeomContainer *tower_geomEM = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_CEMC");
0255 if(!towersIH3 || !towersOH3 || !towersEM3){
0256 std::cout
0257 <<"HCalJetPhiShift::process_event - Error cannot find raw tower node "
0258
0259 << std::endl;
0260 exit(-1);
0261 }
0262
0263 if(!tower_geomIH || !tower_geomOH || !tower_geomEM){
0264 std::cout
0265 <<"HCalJetPhiShift::process_event - Error cannot find raw tower geometry "
0266
0267 << std::endl;
0268 exit(-1);
0269 }
0270
0271 TowerInfo *tower_in, *tower_out;
0272
0273 const int n_channels_IH = (int) towersIH3->size();
0274
0275
0276 for (int i_chan=0; i_chan<n_channels_IH; ++i_chan) {
0277 tower_in = towersIH3->get_tower_at_channel(i_chan);
0278 if (tower_in->get_energy()>0.) {
0279 ++m_nTow_in;
0280
0281 unsigned int calokey = towersIH3->encode_key(i_chan);
0282 int ieta = towersIH3->getTowerEtaBin(calokey);
0283 int iphi = towersIH3->getTowerPhiBin(calokey);
0284
0285 const RawTowerDefs::keytype key = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::HCALIN, ieta, iphi);
0286 float tower_phi = tower_geomIH->get_tower_geometry(key)->get_phi();
0287 float tower_eta = tower_geomIH->get_tower_geometry(key)->get_eta();
0288
0289 m_eta_in.push_back(tower_eta);
0290 m_phi_in.push_back(tower_phi);
0291 m_e_in.push_back(tower_in->get_energy());
0292 m_ieta_in.push_back(ieta);
0293 m_iphi_in.push_back(iphi);
0294 }
0295
0296
0297 tower_out = towersOH3->get_tower_at_channel(i_chan);
0298 if (tower_out->get_energy()>0.) {
0299 ++m_nTow_out;
0300
0301 unsigned int calokey = towersOH3->encode_key(i_chan);
0302 int ieta = towersOH3->getTowerEtaBin(calokey);
0303 int iphi = towersOH3->getTowerPhiBin(calokey);
0304
0305 const RawTowerDefs::keytype key = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::HCALOUT, ieta, iphi);
0306 float tower_phi = tower_geomOH->get_tower_geometry(key)->get_phi();
0307 float tower_eta = tower_geomOH->get_tower_geometry(key)->get_eta();
0308
0309 m_eta_out.push_back(tower_eta);
0310 m_phi_out.push_back(tower_phi);
0311 m_e_out.push_back(tower_out->get_energy());
0312 m_ieta_out.push_back(ieta);
0313 m_iphi_out.push_back(iphi);
0314 }
0315
0316 }
0317
0318 const int n_channels_EM = (int) towersEM3->size();
0319 for (int i_chan=0; i_chan<n_channels_EM; ++i_chan) {
0320
0321
0322 TowerInfo *tower_emc = towersEM3->get_tower_at_channel(i_chan);
0323 if (tower_emc->get_energy()>0.) {
0324 ++m_nTow_emc;
0325
0326 unsigned int calokey = towersEM3->encode_key(i_chan);
0327 int ieta = towersEM3->getTowerEtaBin(calokey);
0328 int iphi = towersEM3->getTowerPhiBin(calokey);
0329
0330 const RawTowerDefs::keytype key = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::CEMC, ieta, iphi);
0331 RawTowerGeom *tower_geom = tower_geomEM->get_tower_geometry(key);
0332 float tower_phi = tower_geom->get_phi();
0333 float tower_eta = tower_geom->get_eta();
0334
0335
0336
0337
0338
0339
0340
0341
0342
0343 m_eta_emc.push_back(tower_eta);
0344 m_phi_emc.push_back(tower_phi);
0345 m_e_emc.push_back(tower_emc->get_energy());
0346 m_ieta_emc.push_back(ieta);
0347 m_iphi_emc.push_back(iphi);
0348 }
0349
0350 }
0351
0352
0353 m_T->Fill();
0354
0355 return Fun4AllReturnCodes::EVENT_OK;
0356 }
