File indexing completed on 2025-08-05 08:13:14
0001 #include "PrintTowers.h"
0002 #include "fastjet/AreaDefinition.hh"
0003 #include "fastjet/ClusterSequenceArea.hh"
0004 #include "fastjet/Selector.hh"
0005 #include "fastjet/tools/BackgroundEstimatorBase.hh"
0006 #include "fastjet/tools/JetMedianBackgroundEstimator.hh"
0007
0008
0009 #include <TFile.h>
0010 #include <TH1F.h>
0011 #include <TH2F.h>
0012 #include <TRandom3.h>
0013 #include <TString.h>
0014 #include <TTree.h>
0015 #include <TH2D.h>
0016 #include <TVector3.h>
0017 #include <algorithm>
0018 #include <cassert>
0019 #include <centrality/CentralityInfo.h>
0020 #include <cmath>
0021 #include <fastjet/ClusterSequence.hh>
0022 #include <fastjet/JetDefinition.hh>
0023 #include <fastjet/PseudoJet.hh>
0024 #include <fun4all/Fun4AllReturnCodes.h>
0025 #include <fun4all/Fun4AllServer.h>
0026 #include <fun4all/PHTFileServer.h>
0027 #include <g4eval/JetEvalStack.h>
0028 #include <jetbase/Jet.h>
0029 #include <jetbase/JetInput.h>
0030 #include <jetbase/TowerJetInput.h>
0031 #include <jetbase/JetMapv1.h>
0032 #include <iostream>
0033 #include <limits>
0034 #include <phool/PHCompositeNode.h>
0035 #include <phool/getClass.h>
0036 #include <stdexcept>
0037 #include <trackbase_historic/SvtxTrackMap.h>
0038 #include <globalvertex/GlobalVertexMap.h>
0039
0040
0041
0042 PrintTowers::PrintTowers
0043 ( std::vector<std::string> jetnames
0044 , std::vector<std::pair<std::string,std::string>> matchlist
0045 , std::vector<std::pair<Jet::SRC,std::string>> src
0046 , const float min_jet_pt
0047 , const float min_tow_pt
0048 )
0049 : SubsysReco("PrintTowers")
0050 , m_jet_names { jetnames }
0051 , m_matchlist { matchlist }
0052 , m_SRC { src }
0053 , m_min_jet_pt { min_jet_pt }
0054 , m_min_tow_pt { min_tow_pt }
0055 { }
0056
0057 PrintTowers::~PrintTowers()
0058 {
0059 }
0060
0061 int PrintTowers::Init(PHCompositeNode* topNode)
0062 { return Fun4AllReturnCodes::EVENT_OK; }
0063
0064 int PrintTowers::End(PHCompositeNode* topNode)
0065 { return Fun4AllReturnCodes::EVENT_OK; }
0066
0067 int PrintTowers::InitRun(PHCompositeNode* topNode)
0068 {
0069 topNode->print();
0070 return Fun4AllReturnCodes::EVENT_OK;
0071 }
0072
0073 int PrintTowers::process_event(PHCompositeNode* topNode)
0074 {
0075
0076 CentralityInfo* cent_node = findNode::getClass<CentralityInfo>(topNode, "CentralityInfo");
0077 float cent = cent_node->get_centile(CentralityInfo::PROP::epd_NS);
0078 if (cent>m_max_cent) return Fun4AllReturnCodes::EVENT_OK;
0079
0080 std::cout << std::endl << " NEW EVENT " << nevnt++ << " cent("<<cent<<")" <<std::endl;
0081
0082
0083
0084 for (auto& name : m_jet_names) {
0085 std::cout << std::endl << " --- JetMap named " << name << " --- " << std::endl;
0086 JetMap* jetmap = findNode::getClass<JetMap>(topNode, name);
0087 if (jetmap == nullptr )
0088 {
0089 std::cout << " -> is NULL (not present) on the node tree." << std::endl;
0090 continue;
0091 }
0092 auto jets = jetmap->vec();
0093 std::sort(jets.begin(), jets.end(), [](Jet* a, Jet* b){return a->get_pt() > b->get_pt();});
0094 std::cout << Form(" -> Has %i jets; printing up to %i jets w/pT>%5.2f", ((int)jets.size()), nmax_jetprint, m_min_jet_pt) << std::endl;
0095 int i {0};
0096 for (auto jet : jets) {
0097 if (jet->get_pt() < m_min_jet_pt) break;
0098 if (i==nmax_jetprint) break;
0099 std::cout << Form("k[%i] pt:phi:eta[%5.2f:%5.2f:%5.2f]", i++, jet->get_pt(), jet->get_phi(), jet->get_eta()) << std::endl;
0100 }
0101 }
0102
0103 for (auto match : m_matchlist) {
0104 JetMap* jetmapA = findNode::getClass<JetMap>(topNode, match.first);
0105 if (jetmapA == nullptr )
0106 {
0107 std::cout << " Can't jet map " << match.first << std::endl;
0108 continue;
0109 }
0110 auto jetsA = jetmapA->vec();
0111 std::sort(jetsA.begin(), jetsA.end(), [](Jet* a, Jet* b){return a->get_pt() > b->get_pt();});
0112
0113 JetMap* jetmapB = findNode::getClass<JetMap>(topNode, match.second);
0114 if (jetmapB == nullptr )
0115 {
0116 std::cout << " Can't find jet map " << match.second << std::endl;
0117 continue;
0118 }
0119 auto jetsB = jetmapB->vec();
0120 std::sort(jetsB.begin(), jetsB.end(), [](Jet* a, Jet* b){return a->get_pt() > b->get_pt();});
0121
0122 std::cout << std::endl << " >>> Matching (up to " << m_pt_min_match << "GeV) " << match.first << " with " << match.second << std::endl;
0123
0124 std::vector<bool> B_is_matched (jetsB.size(),false);
0125 for (auto A : jetsA) {
0126 if (A->get_pt() < m_pt_min_match) break;
0127 bool found_A = false;
0128 for (unsigned int iB = 0; iB<jetsB.size(); ++iB) {
0129 if (B_is_matched[iB]) continue;
0130 auto B = jetsB[iB];
0131 float dphi = fabs(A->get_phi()-B->get_phi());
0132 while (dphi>M_PI) dphi = fabs(dphi - 2*M_PI);
0133 const float deta = A->get_eta() - B->get_eta();
0134 const float R2_comp = deta*deta + dphi*dphi;
0135 if (R2_comp<=0.16) {
0136 std::cout << Form("Matched pT->pT; phi->phi; eta->eta [ %6.2f->%6.2f, %6.2f->%6.2f, %6.2f->%6.2f ]",
0137 A->get_pt(), B->get_pt(), A->get_phi(), B->get_phi(), A->get_eta(), B->get_eta()) << std::endl;
0138 found_A = true;
0139 B_is_matched[iB] = true;
0140 break;
0141 }
0142 }
0143 if (!found_A) {
0144 std::cout << Form("NoMatch pT->??; phi->???? eta->????[ %6.2f->%6s, %6.2f->%6s, %6.2f->%6s ]",
0145 A->get_pt(), "nf", A->get_phi(), "nf", A->get_eta(), "nf") << std::endl;
0146 }
0147 }
0148 }
0149
0150 for (auto src : m_SRC) {
0151 std::cout << std::endl << " --- Printing Tower Input: " << src.second << std::endl;
0152 auto tow_inp = new TowerJetInput(src.first);
0153 auto tows = tow_inp->get_input(topNode);
0154 std::sort(tows.begin(), tows.end(), [](Jet* a, Jet* b){return a->get_pt() > b->get_pt();});
0155 std::cout << Form(" -> Has %i towers; printing all w/pT>%5.2f", ((int)tows.size()), m_min_tow_pt) << std::endl;
0156 int i {0};
0157 for (auto tow : tows) {
0158 if (tow->get_pt() < m_min_tow_pt) break;
0159 std::cout << Form("k[%i] pt:phi:eta[%5.2f:%5.2f:%5.2f]", i++, tow->get_pt(), tow->get_phi(), tow->get_eta()) << std::endl;
0160 }
0161 for (auto tow : tows) delete tow;
0162 }
0163 return Fun4AllReturnCodes::EVENT_OK;
0164
0165 }
0166
