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