File indexing completed on 2025-08-05 08:13:28
0001 #include "TMath.h"
0002 #include "TLorentzVector.h"
0003 #include "MixedJetAnalysis.h"
0004 #include <eventplaneinfo/Eventplaneinfo.h>
0005 #include <eventplaneinfo/EventplaneinfoMap.h>
0006 #include <centrality/CentralityInfov2.h>
0007 #include <calotrigger/MinimumBiasInfov1.h>
0008
0009 #include <TMath.h>
0010 #include <jetbackground/TowerBackground.h>
0011
0012 #include <globalvertex/GlobalVertex.h>
0013 #include <globalvertex/GlobalVertexMap.h>
0014 #include <vector>
0015 #include <mbd/MbdPmtContainer.h>
0016 #include <mbd/MbdPmtHit.h>
0017 #include <fun4all/Fun4AllReturnCodes.h>
0018 #include <phool/PHCompositeNode.h>
0019 #include <phool/PHIODataNode.h>
0020 #include <phool/PHNode.h>
0021 #include <phool/PHNodeIterator.h>
0022 #include <phool/PHObject.h>
0023 #include <phool/getClass.h>
0024
0025
0026 #include <iostream>
0027
0028 #include <map>
0029
0030
0031 MixedJetAnalysis::MixedJetAnalysis(const std::string &name, const std::string &outfilename):
0032 SubsysReco(name)
0033
0034 {
0035 _foutname = outfilename;
0036 }
0037
0038
0039 MixedJetAnalysis::~MixedJetAnalysis()
0040 {
0041
0042 }
0043
0044
0045 int MixedJetAnalysis::Init(PHCompositeNode * )
0046 {
0047
0048
0049
0050 if (Verbosity()) std::cout << __FUNCTION__ << __LINE__<<std::endl;
0051 _f = new TFile( _foutname.c_str(), "RECREATE");
0052
0053 std::cout << " making a file = " << _foutname.c_str() << " , _f = " << _f << std::endl;
0054 _tree_end = new TTree("ttree_end","run stats");
0055 _tree_end->Branch("runnumber", &b_runnumber);
0056 _tree_end->Branch("segment", &b_segment);
0057 _tree_end->Branch("scaled_scalers", &b_scaled_scalers, "scaled_scalers[64]/l");
0058 _tree_end->Branch("live_scalers", &b_live_scalers, "live_scalers[64]/l");
0059 _tree_end->Branch("raw_scalers", &b_raw_scalers, "raw_scalers[64]/l");
0060 _tree_end->Branch("scaled_last", &last_scaled, "scaled_last[64]/l");
0061 _tree_end->Branch("live_last", &last_live, "live_last[64]/l");
0062 _tree_end->Branch("raw_last", &last_raw, "raw_last[64]/l");
0063 _tree_end->Branch("scaled_first", &first_scaled, "scaled_first[64]/l");
0064 _tree_end->Branch("live_first", &first_live, "live_first[64]/l");
0065 _tree_end->Branch("raw_first", &first_raw, "raw_first[64]/l");
0066
0067 _tree = new TTree("ttree","a persevering date tree");
0068
0069 _tree->Branch("gl1_scaled", &b_gl1_scaled, "gl1_scaled/l");
0070 _tree->Branch("gl1_live", &b_gl1_live, "gl1_live/l");
0071 _tree->Branch("gl1_min_bias", &b_gl1_raw, "gl1_min_bias/l");
0072
0073 _i_event = 0;
0074
0075 _tree->Branch("centrality",&b_centrality);
0076 _tree->Branch("psi",&b_psi);
0077 _tree->Branch("v2",&b_v2);
0078 _tree->Branch("flow_fail",&b_flow_fail);
0079 _tree->Branch("minbias",&b_ismb);
0080
0081 _tree->Branch("mbtd_npt", &b_mbd_npmt);
0082 _tree->Branch("mbd_time", &b_mbd_time);
0083 _tree->Branch("mbd_charge", &b_mbd_charge);
0084 _tree->Branch("mbd_side", &b_mbd_side);
0085 _tree->Branch("mbd_ipmt", &b_mbd_ipmt);
0086
0087
0088 _tree->Branch("mbd_vertex_z", &b_vertex_z, "mbd_vertex_z/F");
0089 _tree->Branch("time_zero", &b_time_zero, "time_zero/F");
0090
0091 for (auto cone_size : m_cone_sizes)
0092 {
0093 int coneindex = cone_size.first;
0094 int cone = cone_size.second.first;
0095 bool bkg = cone_size.second.second;
0096 std::cout << " Cone size : " << cone << std::endl;
0097 b_njet.push_back(0);
0098 b_jet_pt.push_back(new std::vector<float>());
0099 b_jet_eta.push_back(new std::vector<float>());
0100 b_jet_phi.push_back(new std::vector<float>());
0101
0102 b_njet_mix.push_back(0);
0103 b_jet_mix_pt.push_back(new std::vector<float>());
0104 b_jet_mix_eta.push_back(new std::vector<float>());
0105 b_jet_mix_phi.push_back(new std::vector<float>());
0106
0107 _tree->Branch(Form("njet_%d%s", cone, (bkg?"_sub":"")), &b_njet.at(coneindex), Form("njet_%d%s/I", cone,(bkg?"_sub":"")));
0108 _tree->Branch(Form("jet_pt_%d%s", cone, (bkg?"_sub":"")), b_jet_pt.at(coneindex));
0109 _tree->Branch(Form("jet_eta_%d%s", cone, (bkg?"_sub":"")), b_jet_eta.at(coneindex));
0110 _tree->Branch(Form("jet_phi_%d%s", cone, (bkg?"_sub":"")), b_jet_phi.at(coneindex));
0111
0112 _tree->Branch(Form("njet_mix_%d%s", cone, (bkg?"_sub":"")), &b_njet_mix.at(coneindex), Form("njet_mix_%d%s/I", cone,(bkg?"_sub":"")));
0113 _tree->Branch(Form("jet_mix_pt_%d%s", cone, (bkg?"_sub":"")), b_jet_mix_pt.at(coneindex));
0114 _tree->Branch(Form("jet_mix_eta_%d%s", cone, (bkg?"_sub":"")), b_jet_mix_eta.at(coneindex));
0115 _tree->Branch(Form("jet_mix_phi_%d%s", cone, (bkg?"_sub":"")), b_jet_mix_phi.at(coneindex));
0116
0117 std::map<unsigned int, std::vector<std::vector<struct mixedjet>>*> empty_map;
0118 unsigned int ncentbins = (m_use_psi?20:100);
0119 for (unsigned int i = 0; i < ncentbins; i++)
0120 {
0121 for (unsigned int j = 0; j < 12; j++)
0122 {
0123 if (m_use_psi)
0124 {
0125 for (unsigned int k = 0; k < 8; k++)
0126 {
0127 unsigned int total_bin = (i & 0x1f) + ((j & 0xf) << 0x5U) + ((k & 0x7) << 0x9U);
0128
0129 empty_map[total_bin] = new std::vector<std::vector<struct mixedjet>>;;
0130 }
0131 }
0132 else
0133 {
0134 unsigned int total_bin = (i & 0xff) + ((j & 0xf) << 0x8U);
0135
0136 empty_map[total_bin] = new std::vector<std::vector<struct mixedjet>>;;
0137
0138 }
0139 }
0140 }
0141 m_mixed_jet_bank.push_back(empty_map);
0142 }
0143 std::cout << "Done initing the treemaker"<<std::endl;
0144 return Fun4AllReturnCodes::EVENT_OK;
0145 }
0146
0147
0148
0149 void MixedJetAnalysis::reset_tree_vars()
0150 {
0151 if (Verbosity()) std::cout << __FUNCTION__ << __LINE__<<std::endl;
0152
0153 b_ismb = 0;
0154 b_centrality = -99;
0155 b_mbd_npmt = 0;
0156 b_mbd_charge.clear();
0157 b_mbd_time.clear();
0158 b_mbd_ipmt.clear();
0159 b_mbd_side.clear();
0160
0161 for (int i = 0; i < m_n_cone_sizes; i++)
0162 {
0163 b_njet.at(i) = 0;
0164 b_jet_pt.at(i)->clear();
0165 b_jet_eta.at(i)->clear();
0166 b_jet_phi.at(i)->clear();
0167
0168 b_njet_mix.at(i) = 0;
0169 b_jet_mix_pt.at(i)->clear();
0170 b_jet_mix_eta.at(i)->clear();
0171 b_jet_mix_phi.at(i)->clear();
0172 }
0173
0174 b_gl1_scaled = 0x0;
0175 b_gl1_live = 0x0;
0176 b_gl1_raw = 0x0;
0177
0178
0179 return;
0180 }
0181
0182 int MixedJetAnalysis::process_event(PHCompositeNode *topNode)
0183 {
0184
0185 if (Verbosity()) std::cout << __FILE__ << " "<< __LINE__<<" Event : " << _i_event <<std::endl;
0186
0187 reset_tree_vars();
0188
0189 Gl1Packet* gl1p = findNode::getClass<Gl1Packet>(topNode, 14001);
0190 b_gl1_scaled = gl1p->getScaledVector();
0191 b_gl1_live = gl1p->getLiveVector();
0192 b_gl1_raw = gl1p->lValue(10, 1);
0193
0194 if (_i_event == 0)
0195 {
0196 for (int i = 0; i < 64; i++)
0197 {
0198 first_raw[i] = gl1p->lValue(i, 0);
0199 first_live[i] = gl1p->lValue(i, 1);
0200 first_scaled[i] = gl1p->lValue(i, 2);
0201 }
0202 }
0203 else
0204 {
0205 for (int i = 0; i < 64; i++)
0206 {
0207 last_raw[i] = gl1p->lValue(i, 0);
0208 last_live[i] = gl1p->lValue(i, 1);
0209 last_scaled[i] = gl1p->lValue(i, 2);
0210 }
0211 }
0212
0213 _i_event++;
0214
0215
0216
0217
0218
0219
0220 m_minimumbiasinfo = findNode::getClass<MinimumBiasInfo>(topNode, "MinimumBiasInfo");
0221 m_towback = findNode::getClass<TowerBackground>(topNode, "TowerInfoBackground_Sub2");
0222 m_central = findNode::getClass<CentralityInfov2>(topNode, "CentralityInfo");
0223 m_evpmap = findNode::getClass<EventplaneinfoMap>(topNode, "EventplaneinfoMap");
0224
0225 b_centrality = 0;
0226 b_ismb = 0;
0227 b_psi = 0;
0228 b_v2 = -99;
0229 if (!m_evpmap && m_use_psi)
0230 {
0231 exit(-1);
0232 }
0233 if (m_use_psi)
0234 {
0235 if(!(m_evpmap->empty()))
0236 {
0237
0238 auto EPDNS = m_evpmap->get(EventplaneinfoMap::sEPDNS);
0239 b_psi = EPDNS->get_shifted_psi(2);
0240 }
0241 }
0242 if (m_central)
0243 b_centrality = (m_central->has_centile(CentralityInfo::PROP::mbd_NS)?m_central->get_centrality_bin(CentralityInfo::PROP::mbd_NS) : -99);
0244 if (m_minimumbiasinfo)
0245 b_ismb = (m_minimumbiasinfo->isAuAuMinimumBias()? 1 : 0);
0246 if (m_towback)
0247 {
0248 b_v2 = m_towback->get_v2();
0249 b_flow_fail = m_towback->get_flow_failure_flag();
0250 }
0251 GlobalVertexMap* vertexmap = findNode::getClass<GlobalVertexMap>(topNode, "GlobalVertexMap");
0252
0253
0254 float vtx_z = 0;
0255 b_vertex_z = -999;
0256 if (vertexmap && !vertexmap->empty())
0257 {
0258 GlobalVertex* vtx = vertexmap->begin()->second;
0259 if (vtx)
0260 {
0261 vtx_z = vtx->get_z();
0262 b_vertex_z = vtx_z;
0263 b_time_zero = vtx->get_t();
0264 }
0265 }
0266 if (Verbosity()) std::cout << "Mbd Vertex = " << b_vertex_z << std::endl;
0267
0268 if (!b_ismb) return Fun4AllReturnCodes::EVENT_OK;
0269
0270 if (Verbosity()) std::cout << "MINBIAS event"<< std::endl;
0271
0272 unsigned int centrality_bin = b_centrality/5;
0273
0274 unsigned int vertex_bin = 0;
0275
0276 for (unsigned int i = 0; i < m_nbin_vertex; i++)
0277 {
0278 if (b_vertex_z < m_vertex_bin_edges[i+1])
0279 {
0280 vertex_bin = i;
0281 break;
0282 }
0283 }
0284 unsigned int psi_bin = 0;
0285 if (m_use_psi)
0286 {
0287 for (unsigned int i = 0; i < m_nbin_psi; i++)
0288 {
0289 if (b_psi < (-1*TMath::Pi()/2. + (i+1)+TMath::Pi()/8.))
0290 {
0291 psi_bin = i;
0292 break;
0293 }
0294 }
0295 m_total_bin = (centrality_bin & 0x1f) + ((vertex_bin & 0xf) << 0x5U) + ((psi_bin & 0x7) << 9U);
0296 }
0297 else
0298 {
0299 m_total_bin = (centrality_bin & 0xff) + ((vertex_bin & 0xf) << 0x8U);
0300 }
0301
0302
0303 if (Verbosity()) std::cout << " total bin : " << m_total_bin << std::endl;
0304 bool keep_candidate = false;
0305 for (int i = 0; i < m_n_cone_sizes; i++)
0306 {
0307 bool jet_candidate = process_jets(i, topNode);
0308
0309 if (jet_candidate)
0310 {
0311 if (Verbosity()) std::cout << "Found good Jet" << std::endl;
0312 keep_candidate = true;
0313 get_mixed_events(i);
0314 }
0315 if (Verbosity()) std::cout << "Now storing jets vector" << std::endl;
0316 store_mixed_events(i);
0317
0318 if (Verbosity())
0319 {
0320 show_queue(i);
0321 }
0322 }
0323
0324 if (!keep_candidate) return Fun4AllReturnCodes::EVENT_OK;
0325
0326 MbdPmtContainer *pmts_mbd = findNode::getClass<MbdPmtContainer>(topNode, "MbdPmtContainer");
0327
0328 if (pmts_mbd)
0329 {
0330
0331 for (int i = 0 ; i < 128; i++)
0332 {
0333 MbdPmtHit *tmp_pmt = pmts_mbd->get_pmt(i);
0334
0335 float charge = tmp_pmt->get_q();
0336 float time = tmp_pmt->get_time();
0337 if (fabs(time) < 25. && charge > 0.4)
0338 {
0339 b_mbd_charge.push_back(charge);
0340 b_mbd_time.push_back(time);
0341 b_mbd_ipmt.push_back(i);
0342 b_mbd_side.push_back(i/64);
0343 b_mbd_npmt++;
0344 }
0345 }
0346 }
0347
0348 _tree->Fill();
0349
0350 return Fun4AllReturnCodes::EVENT_OK;
0351 }
0352
0353 int MixedJetAnalysis::process_jets(int cone_index, PHCompositeNode* topNode)
0354 {
0355
0356 bool isbkg = m_cone_sizes[cone_index].second;
0357 int cone_size = m_cone_sizes[cone_index].first;
0358
0359 std::string recoJetName = Form("AntiKt_Tower_r%02d", cone_size);
0360
0361 bool jet_candidate = false;
0362
0363 if (isbkg)
0364 {
0365 recoJetName = Form("AntiKt_Tower_r%02d_Sub1", cone_size);
0366 }
0367
0368 JetContainer *jets = findNode::getClass<JetContainer>(topNode, recoJetName);
0369 if (jets)
0370 {
0371
0372 for (auto jet : *jets)
0373 {
0374
0375 if (jet->get_pt() < pt_cut) continue;
0376
0377 if (jet->get_pt() >= signal_pt_cut) jet_candidate = true;
0378
0379 b_njet.at(cone_index)++;
0380 b_jet_pt.at(cone_index)->push_back(jet->get_pt());
0381 b_jet_eta.at(cone_index)->push_back(jet->get_eta());
0382 b_jet_phi.at(cone_index)->push_back(jet->get_phi());
0383
0384 }
0385 }
0386
0387 return jet_candidate;
0388
0389 }
0390
0391 int MixedJetAnalysis::store_mixed_events(int cone_index)
0392 {
0393
0394 auto mixed_jet_vector = m_mixed_jet_bank.at(cone_index)[m_total_bin];
0395 if (Verbosity()) std::cout << "Storing Mixed events" << std::endl;
0396 if (!mixed_jet_vector) return 1;
0397 if (Verbosity()) std::cout << "entering Mixed event" << std::endl;
0398
0399 int njet = b_njet.at(cone_index);
0400
0401 if (mixed_jet_vector->size() > m_max_length_buffer) return 0;
0402 std::vector<struct mixedjet> moxies = {};
0403 for (int ij = 0 ; ij < njet ; ij++)
0404 {
0405 struct mixedjet moxy;
0406 moxy.pt = b_jet_pt.at(cone_index)->at(ij);
0407 moxy.eta = b_jet_eta.at(cone_index)->at(ij);
0408 moxy.phi = b_jet_phi.at(cone_index)->at(ij);
0409 moxies.push_back(moxy);
0410 }
0411
0412 mixed_jet_vector->push_back(moxies);
0413
0414 return 0;
0415
0416 }
0417
0418 int MixedJetAnalysis::get_mixed_events(int cone_index)
0419 {
0420
0421 auto mixed_jet_vector = m_mixed_jet_bank.at(cone_index)[m_total_bin];
0422 if (Verbosity()) std::cout << "Getting Mixed events" << std::endl;
0423 if (!mixed_jet_vector) return 1;
0424 if (Verbosity()) std::cout << "Got Mixed event" << std::endl;
0425 int njet = mixed_jet_vector->size();
0426 if (Verbosity()) std::cout << "Found " << njet << " Jet vectors" << std::endl;
0427
0428 if (njet == 0)
0429 {
0430 b_njet_mix.at(cone_index) = -1;
0431 return 1;
0432 }
0433 auto moxies = mixed_jet_vector->at(0);
0434
0435 b_njet_mix.at(cone_index) = 0;
0436 for (auto moxy : moxies)
0437 {
0438 if (Verbosity()) std::cout << " getting mixed jet: " << moxy.pt << " / " << moxy.eta << " / " << moxy.phi << std::endl;
0439 b_njet_mix.at(cone_index)++;
0440 b_jet_mix_pt.at(cone_index)->push_back(moxy.pt);
0441 b_jet_mix_eta.at(cone_index)->push_back(moxy.eta);
0442 b_jet_mix_phi.at(cone_index)->push_back(moxy.phi);
0443 }
0444 if (Verbosity()) std::cout << "Now erasing vector" << std::endl;
0445 mixed_jet_vector->erase(mixed_jet_vector->begin());
0446
0447 return 0;
0448
0449 }
0450
0451 int MixedJetAnalysis::show_queue(int cone_index)
0452 {
0453 for (unsigned int i = 0; i < 100; i++)
0454 {
0455 std::cout << i << "\t | ";
0456 for (unsigned int j = 0; j < 1; j++)
0457 {
0458 unsigned int total_bin = (i & 0xff) + ((j & 0xff) << 0x8U);
0459 int number_in_queue = m_mixed_jet_bank.at(cone_index)[total_bin]->size();
0460 std::cout << number_in_queue << " \t";
0461 }
0462 std::cout << " | " << std::endl;
0463 }
0464 return 0;
0465 }
0466
0467 int MixedJetAnalysis::End(PHCompositeNode *)
0468 {
0469 if (Verbosity() > 0)
0470 {
0471 std::cout << "MixedJetAnalysis::End(PHCompositeNode *topNode) This is the End..." << std::endl;
0472 }
0473 std::cout<<"Total events: "<<_i_event<<std::endl;
0474
0475 for (int i = 0; i < 64; i++)
0476 {
0477 b_scaled_scalers[i] = (last_scaled[i] - first_scaled[i]);
0478 b_live_scalers[i] = (last_live[i] - first_live[i]);
0479 b_raw_scalers[i] = (last_raw[i] - first_raw[i]);
0480 }
0481
0482 _tree_end->Fill();
0483
0484 _f->Write();
0485 _f->Close();
0486
0487 return Fun4AllReturnCodes::EVENT_OK;
0488 }
0489
