File indexing completed on 2025-08-06 08:17:44
0001 #include "JetBackgroundCut.h"
0002
0003 #include <calobase/RawTowerDefs.h> // for CalorimeterId, encode_to...
0004 #include <calobase/RawTowerGeom.h>
0005 #include <calobase/RawTowerGeomContainer.h>
0006 #include <calobase/TowerInfo.h> // for TowerInfo
0007 #include <calobase/TowerInfoContainer.h>
0008
0009 #include <globalvertex/GlobalVertexMap.h>
0010 #include <globalvertex/Vertex.h> // for Vertex
0011
0012 #include <jetbase/Jet.h>
0013 #include <jetbase/JetContainer.h>
0014
0015 #include <phparameter/PHParameters.h>
0016
0017 #include <fun4all/Fun4AllReturnCodes.h>
0018
0019 #include <phool/PHCompositeNode.h>
0020 #include <phool/getClass.h>
0021
0022 #include <cmath>
0023 #include <iostream> // for basic_ostream, operator<<
0024 #include <map> // for _Rb_tree_iterator, opera...
0025 #include <utility> // for pair
0026 #include <vector> // for vector
0027
0028
0029 JetBackgroundCut::JetBackgroundCut(const std::string &jetNodeName, const std::string &name, const int debug, const bool doAbort, GlobalVertex::VTXTYPE vtxtype, int sysvar)
0030 : SubsysReco(name)
0031 , _doAbort(doAbort)
0032 , _name(name)
0033 , _debug(debug)
0034 , _jetNodeName(jetNodeName)
0035 , _vtxtype(vtxtype)
0036 , _sysvar(sysvar)
0037 , _cutParams(name)
0038 {
0039 SetDefaultParams();
0040 }
0041
0042
0043 int JetBackgroundCut::Init(PHCompositeNode *topNode)
0044 {
0045 CreateNodeTree(topNode);
0046
0047 return Fun4AllReturnCodes::EVENT_OK;
0048 }
0049
0050 void JetBackgroundCut::CreateNodeTree(PHCompositeNode *topNode)
0051 {
0052 PHNodeIterator iter(topNode);
0053
0054 PHCompositeNode *parNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "PAR"));
0055 if (!parNode)
0056 {
0057 std::cout << "No RUN node found; cannot create PHParameters for storing cut results. Aborting run!";
0058 }
0059
0060 _cutParams.SaveToNodeTree(parNode, "JetCutParams");
0061 }
0062
0063
0064 int JetBackgroundCut::process_event(PHCompositeNode *topNode)
0065 {
0066 TowerInfoContainer *towersEM = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_CEMC_RETOWER");
0067 TowerInfoContainer *towersOH = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_HCALOUT");
0068 JetContainer *jets = findNode::getClass<JetContainer>(topNode, _jetNodeName);
0069 GlobalVertexMap *gvtxmap = findNode::getClass<GlobalVertexMap>(topNode, "GlobalVertexMap");
0070
0071 RawTowerGeomContainer *geom[2];
0072 geom[0] = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALIN");
0073 geom[1] = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALOUT");
0074
0075 float zvtx = NAN;
0076 float maxJetET = 0;
0077 float maxJetPhi = NAN;
0078 float subJetET = 0;
0079 float subJetPhi = NAN;
0080 float frcem = 0;
0081 float frcoh = 0;
0082 float dPhi = NAN;
0083
0084 if (!towersEM || !towersOH || !geom[0] || !geom[1] || !gvtxmap)
0085 {
0086 if (_debug > 0 && !_missingInfoWarningPrinted)
0087 {
0088 std::cout << "Missing critical info; abort event. Further warnings will be suppressed. AddressOf towersEM/towersOH/geomIH/geomOH/gvtxmap : " << towersEM << "/" << towersOH << "/" << geom[0] << "/" << geom[1] << "/" << gvtxmap << std::endl;
0089 }
0090 _missingInfoWarningPrinted = true;
0091 return Fun4AllReturnCodes::ABORTEVENT;
0092 }
0093
0094 if (gvtxmap)
0095 {
0096 GlobalVertex *gvtx = nullptr;
0097 if (gvtxmap->empty())
0098 {
0099 if (_debug > 0)
0100 {
0101 std::cout << "gvtxmap empty - set zvtx to 0 and continue." << std::endl;
0102 }
0103 zvtx = 0;
0104 }
0105 else
0106 {
0107 gvtx = gvtxmap->begin()->second;
0108 }
0109 if (gvtx)
0110 {
0111 auto startIter = gvtx->find_vertexes(_vtxtype);
0112 auto endIter = gvtx->end_vertexes();
0113 for (auto iter = startIter; iter != endIter; ++iter)
0114 {
0115 const auto &[type, vertexVec] = *iter;
0116 if (type != _vtxtype)
0117 {
0118 continue;
0119 }
0120 for (const auto *vertex : vertexVec)
0121 {
0122 if (!vertex)
0123 {
0124 continue;
0125 }
0126 zvtx = vertex->get_z();
0127 }
0128 }
0129 }
0130 else
0131 {
0132 if (_debug > 0)
0133 {
0134 std::cout << "gvtx is NULL! Set zvtx to 0 and continue." << std::endl;
0135 }
0136 zvtx = 0;
0137 }
0138 }
0139 else
0140 {
0141 if (_debug > 0)
0142 {
0143 std::cout << "gvtxmap is NULL! ABORT EVENT!" << std::endl;
0144 }
0145 return Fun4AllReturnCodes::ABORTEVENT;
0146 }
0147
0148 if (std::isnan(zvtx))
0149 {
0150 if (_debug > 0)
0151 {
0152 std::cout << "zvtx is NAN after attempting to grab it. Set to 0 and continue." << std::endl;
0153 }
0154 zvtx = 0;
0155 }
0156
0157 if (_debug > 1)
0158 {
0159 std::cout << "Getting jets: " << std::endl;
0160 }
0161
0162 if (jets)
0163 {
0164 int tocheck = jets->size();
0165 if (_debug > 2)
0166 {
0167 std::cout << "Found " << tocheck << " jets to check..." << std::endl;
0168 }
0169 for (int i = 0; i < tocheck; ++i)
0170 {
0171 float jetET = 0;
0172 float jetPhi = NAN;
0173 float jetEta = NAN;
0174 Jet *jet = jets->get_jet(i);
0175 if (jet)
0176 {
0177 jetEta = jet->get_eta();
0178 jetET = jet->get_e() / std::cosh(jetEta);
0179 jetPhi = jet->get_phi();
0180 }
0181 else
0182 {
0183 continue;
0184 }
0185 if (jetET < 8)
0186 {
0187 continue;
0188 }
0189 if (_debug > 2)
0190 {
0191 std::cout << "found a good jet!" << std::endl;
0192 }
0193 if (jetET > maxJetET)
0194 {
0195 if (maxJetET)
0196 {
0197 subJetET = maxJetET;
0198 subJetPhi = maxJetPhi;
0199 }
0200 maxJetET = jetET;
0201 maxJetPhi = jetPhi;
0202 }
0203 else if (jetET > subJetET)
0204 {
0205 subJetET = jetET;
0206 subJetPhi = jetPhi;
0207 continue;
0208 }
0209 else
0210 {
0211 continue;
0212 }
0213 frcem = 0;
0214 frcoh = 0;
0215 for (auto comp : jet->get_comp_vec())
0216 {
0217 unsigned int channel = comp.second;
0218 TowerInfo *tower;
0219 if (comp.first == 13 || comp.first == 28 || comp.first == 25)
0220 {
0221 tower = towersEM->get_tower_at_channel(channel);
0222 int key = towersEM->encode_key(channel);
0223 const RawTowerDefs::keytype geomkey = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::HCALIN, towersEM->getTowerEtaBin(key), towersEM->getTowerPhiBin(key));
0224 RawTowerGeom *tower_geom = geom[0]->get_tower_geometry(geomkey);
0225 float radius = 93.5;
0226 float ihEta = tower_geom->get_eta();
0227 float emZ = radius / (std::tan(2 * std::atan(std::exp(-ihEta))));
0228 float newz = emZ - zvtx;
0229 float newTheta = std::atan2(radius, newz);
0230 float towerEta = -log(std::tan(0.5 * newTheta));
0231 frcem += tower->get_energy() / std::cosh(towerEta);
0232 }
0233 if (comp.first == 7 || comp.first == 27)
0234 {
0235 tower = towersOH->get_tower_at_channel(channel);
0236 int key = towersOH->encode_key(channel);
0237 const RawTowerDefs::keytype geomkey = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::HCALOUT, towersOH->getTowerEtaBin(key), towersOH->getTowerPhiBin(key));
0238 RawTowerGeom *tower_geom = geom[1]->get_tower_geometry(geomkey);
0239 float radius = tower_geom->get_center_radius();
0240 float newz = tower_geom->get_center_z() - zvtx;
0241 float newTheta = std::atan2(radius, newz);
0242 float towerEta = -log(std::tan(0.5 * newTheta));
0243 frcoh += tower->get_energy() / std::cosh(towerEta);
0244 }
0245 }
0246
0247 frcem /= maxJetET;
0248 frcoh /= maxJetET;
0249 }
0250 }
0251 else
0252 {
0253 if (_debug > 0)
0254 {
0255 std::cout << "No jet node!" << std::endl;
0256 }
0257 return Fun4AllReturnCodes::ABORTEVENT;
0258 }
0259
0260 bool isDijet = false;
0261 if (subJetET > 8)
0262 {
0263 isDijet = true;
0264 dPhi = std::abs(maxJetPhi - subJetPhi);
0265 if (dPhi > M_PI)
0266 {
0267 dPhi = 2 * M_PI - dPhi;
0268 }
0269 }
0270 bool dPhiCut = failsdPhiCut(dPhi, isDijet);
0271
0272 bool failsLoEm = failsLoEmFracETCut(frcem, maxJetET, dPhiCut, isDijet);
0273 bool failsHiEm = failsHiEmFracETCut(frcem, maxJetET, dPhiCut, isDijet);
0274 bool failsIhCut = failsIhFracCut(frcem, frcoh);
0275
0276 bool failsAnyCut = failsLoEm || failsHiEm || failsIhCut;
0277
0278 if (failsAnyCut && _doAbort)
0279 {
0280 return Fun4AllReturnCodes::ABORTEVENT;
0281 }
0282 PHNodeIterator iter(topNode);
0283 PHCompositeNode *parNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "PAR"));
0284 _cutParams.set_int_param("failsLoEmJetCut", failsLoEm);
0285 _cutParams.set_int_param("failsHiEmJetCut", failsHiEm);
0286 _cutParams.set_int_param("failsIhJetCut", failsIhCut);
0287 _cutParams.set_int_param("failsAnyJetCut", failsAnyCut);
0288 _cutParams.set_int_param("isDijet", isDijet);
0289 _cutParams.set_double_param("frcem", frcem);
0290 _cutParams.set_double_param("frcoh", frcoh);
0291 _cutParams.set_double_param("maxJetET", maxJetET);
0292 _cutParams.set_double_param("dPhi", dPhi);
0293 _cutParams.UpdateNodeTree(parNode, "JetCutParams");
0294
0295 return Fun4AllReturnCodes::EVENT_OK;
0296 }
0297
0298 int JetBackgroundCut::ResetEvent(PHCompositeNode * )
0299 {
0300 if (Verbosity() > 0)
0301 {
0302 std::cout << "JetBackgroundCut::ResetEvent(PHCompositeNode *topNode) Resetting internal structures, prepare for next event" << std::endl;
0303 }
0304 return Fun4AllReturnCodes::EVENT_OK;
0305 }
0306
0307
0308 int JetBackgroundCut::End(PHCompositeNode * )
0309 {
0310 if (Verbosity() > 0)
0311 {
0312 std::cout << "JetBackgroundCut::End(PHCompositeNode *topNode) This is the End..." << std::endl;
0313 }
0314 return Fun4AllReturnCodes::EVENT_OK;
0315 }
0316
0317
0318 int JetBackgroundCut::Reset(PHCompositeNode * )
0319 {
0320 if (Verbosity() > 0)
0321 {
0322 std::cout << "JetBackgroundCut::Reset(PHCompositeNode *topNode) being Reset" << std::endl;
0323 }
0324 return Fun4AllReturnCodes::EVENT_OK;
0325 }
0326
0327
0328 void JetBackgroundCut::Print(const std::string &what) const
0329 {
0330 std::cout << "JetBackgroundCut::Print(const std::string &what) const Printing info for " << what << std::endl;
0331 }
