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