File indexing completed on 2025-08-06 08:17:36
0001 #include "CentralityReco.h"
0002
0003 #include "CentralityInfov2.h"
0004
0005 #include <mbd/MbdOut.h>
0006 #include <mbd/MbdPmtContainer.h>
0007 #include <mbd/MbdPmtHit.h>
0008
0009 #include <calotrigger/MinimumBiasInfo.h>
0010
0011 #include <globalvertex/GlobalVertexMap.h>
0012 #include <globalvertex/GlobalVertex.h>
0013
0014 #include <ffamodules/CDBInterface.h>
0015 #include <cdbobjects/CDBTTree.h>
0016
0017 #include <fun4all/Fun4AllReturnCodes.h>
0018
0019 #include <phool/PHCompositeNode.h>
0020 #include <phool/PHIODataNode.h>
0021 #include <phool/PHNode.h>
0022 #include <phool/PHNodeIterator.h>
0023 #include <phool/PHObject.h>
0024 #include <phool/getClass.h>
0025 #include <phool/phool.h>
0026
0027 #include <cmath>
0028 #include <filesystem>
0029 #include <iostream>
0030 #include <string>
0031
0032 CentralityReco::CentralityReco(const std::string &name)
0033 : SubsysReco(name)
0034 {
0035
0036 }
0037
0038 int CentralityReco::InitRun(PHCompositeNode *topNode)
0039 {
0040 CDBInterface *m_cdb = CDBInterface::instance();
0041
0042 std::string centdiv_url = m_cdb->getUrl("Centrality");
0043 if (m_overwrite_divs)
0044 {
0045 centdiv_url = m_overwrite_url_divs;
0046 std::cout << " Overwriting Divs to " << m_overwrite_url_divs << std::endl;
0047 }
0048
0049 if (Download_centralityDivisions(centdiv_url))
0050 {
0051 return Fun4AllReturnCodes::ABORTRUN;
0052 }
0053
0054 std::string centscale_url = m_cdb->getUrl("CentralityScale");
0055 if (m_overwrite_scale)
0056 {
0057 centscale_url = m_overwrite_url_scale;
0058 std::cout << " Overwriting Scale to " << m_overwrite_url_scale << std::endl;
0059 }
0060
0061 if (Download_centralityScale(centscale_url))
0062 {
0063 return Fun4AllReturnCodes::ABORTRUN;
0064 }
0065
0066 std::string vertexscale_url = m_cdb->getUrl("CentralityVertexScale");
0067 if (m_overwrite_vtx)
0068 {
0069 vertexscale_url = m_overwrite_url_vtx;
0070 std::cout << " Overwriting Vtx to " << m_overwrite_url_vtx << std::endl;
0071 }
0072
0073 if (Download_centralityVertexScales(vertexscale_url))
0074 {
0075 return Fun4AllReturnCodes::ABORTRUN;
0076 }
0077
0078 CreateNodes(topNode);
0079 return Fun4AllReturnCodes::EVENT_OK;
0080 }
0081
0082 int CentralityReco::Download_centralityScale(const std::string &dbfile)
0083 {
0084 m_centrality_scale = 1.00;
0085
0086 std::filesystem::path dbase_file = dbfile;
0087 if (dbase_file.extension() == ".root")
0088 {
0089 CDBTTree *cdbttree = new CDBTTree(dbase_file);
0090 cdbttree->LoadCalibrations();
0091 m_centrality_scale = cdbttree->GetDoubleValue(0, "centralityscale");
0092 if (Verbosity())
0093 {
0094 std::cout << "centscale = " << m_centrality_scale << std::endl;
0095 }
0096 delete cdbttree;
0097 }
0098 else
0099 {
0100 std::cout << PHWHERE << ", ERROR, unknown file type, " << dbfile << std::endl;
0101 return Fun4AllReturnCodes::ABORTRUN;
0102 }
0103
0104 return Fun4AllReturnCodes::EVENT_OK;
0105 }
0106
0107 int CentralityReco::Download_centralityDivisions(const std::string &dbfile)
0108 {
0109 m_centrality_map.fill(0);
0110
0111 std::filesystem::path dbase_file = dbfile;
0112
0113 if (dbase_file.extension() == ".root")
0114 {
0115 CDBTTree *cdbttree = new CDBTTree(dbase_file);
0116 cdbttree->LoadCalibrations();
0117 if (Verbosity())
0118 {
0119 cdbttree->Print();
0120 }
0121 for (int idiv = 0; idiv < NDIVS; idiv++)
0122 {
0123 m_centrality_map[idiv] = cdbttree->GetFloatValue(idiv, "centralitydiv");
0124 if (Verbosity())
0125 {
0126 std::cout << "centdiv " << idiv << " : " << m_centrality_map[idiv] << std::endl;
0127 }
0128 }
0129 delete cdbttree;
0130 }
0131 else
0132 {
0133 std::cout << PHWHERE << ", ERROR, unknown file type, " << dbfile << std::endl;
0134 return Fun4AllReturnCodes::ABORTRUN;
0135 }
0136
0137 return Fun4AllReturnCodes::EVENT_OK;
0138 }
0139 int CentralityReco::Download_centralityVertexScales(const std::string &dbfile)
0140 {
0141
0142 std::filesystem::path dbase_file = dbfile;
0143
0144 if (dbase_file.extension() == ".root")
0145 {
0146 CDBTTree *cdbttree = new CDBTTree(dbase_file);
0147 cdbttree->LoadCalibrations();
0148 if (Verbosity())
0149 {
0150 cdbttree->Print();
0151 }
0152
0153 int nvertexbins = cdbttree->GetIntValue(0, "nvertexbins");
0154
0155
0156 for (int iv = 0; iv < nvertexbins; iv++)
0157 {
0158 float scale = cdbttree->GetDoubleValue(iv, "scale");
0159 float lowvertex = cdbttree->GetDoubleValue(iv, "low_vertex");
0160 float highvertex = cdbttree->GetDoubleValue(iv, "high_vertex");
0161 m_vertex_scales.emplace_back(std::make_pair(lowvertex, highvertex), scale);
0162 }
0163
0164 delete cdbttree;
0165 }
0166 else
0167 {
0168 std::cout << PHWHERE << ", ERROR, unknown file type, " << dbfile << std::endl;
0169 return Fun4AllReturnCodes::ABORTRUN;
0170 }
0171
0172 return Fun4AllReturnCodes::EVENT_OK;
0173 }
0174
0175 int CentralityReco::ResetEvent(PHCompositeNode * )
0176 {
0177 m_mbd_total_charge = 0;
0178
0179 return Fun4AllReturnCodes::EVENT_OK;
0180 }
0181
0182 int CentralityReco::FillVars()
0183 {
0184 if (Verbosity() > 1)
0185 {
0186 std::cout << __FILE__ << " :: " << __FUNCTION__ << std::endl;
0187 }
0188
0189
0190 float scale_factor = getVertexScale();
0191
0192 if (Verbosity())
0193 {
0194 std::cout << scale_factor << "*" << m_centrality_scale << std::endl;
0195 }
0196
0197 for (int i = 0; i < 128; i++)
0198 {
0199
0200 m_mbd_hit = m_mbd_container->get_pmt(i);
0201
0202 if ((m_mbd_hit->get_q()) < mbd_charge_cut)
0203 {
0204 continue;
0205 }
0206 if (fabs(m_mbd_hit->get_time()) > mbd_time_cut)
0207 {
0208 continue;
0209 }
0210 m_mbd_total_charge += m_mbd_hit->get_q()*scale_factor*m_centrality_scale;
0211 }
0212
0213
0214 return Fun4AllReturnCodes::EVENT_OK;
0215 }
0216
0217 int CentralityReco::FillCentralityInfo()
0218 {
0219
0220
0221 if (Verbosity() > 1)
0222 {
0223 std::cout << __FILE__ << " :: " << __FUNCTION__ << std::endl;
0224 }
0225
0226 int binvalue = std::numeric_limits<int>::quiet_NaN();
0227 float value = std::numeric_limits<float>::quiet_NaN();
0228 for (int i = 0; i < NDIVS; i++)
0229 {
0230 if (m_centrality_map[i] < m_mbd_total_charge)
0231 {
0232 binvalue = i + 1;
0233 value = static_cast<float>(i + 1)/static_cast<float>(NDIVS);
0234 break;
0235 }
0236 }
0237 if (Verbosity())
0238 {
0239 std::cout << " Centile : " << value << std::endl;
0240 std::cout << " Charge : " << m_mbd_total_charge << std::endl;
0241 }
0242
0243 m_central->set_centile(CentralityInfo::PROP::mbd_NS, value);
0244 m_central->set_centrality_bin(CentralityInfo::PROP::mbd_NS, binvalue);
0245
0246 return Fun4AllReturnCodes::EVENT_OK;
0247 }
0248
0249 int CentralityReco::process_event(PHCompositeNode *topNode)
0250 {
0251 if (Verbosity())
0252 {
0253 std::cout << "------------CentralityReco-------------" << std::endl;
0254 }
0255
0256
0257 if (GetNodes(topNode))
0258 {
0259 return Fun4AllReturnCodes::ABORTRUN;
0260 }
0261
0262 if (!m_mb_info->isAuAuMinimumBias())
0263 {
0264 return Fun4AllReturnCodes::EVENT_OK;
0265 }
0266
0267
0268
0269 if (FillVars())
0270 {
0271 return Fun4AllReturnCodes::ABORTEVENT;
0272 }
0273
0274 if (FillCentralityInfo())
0275 {
0276 return Fun4AllReturnCodes::ABORTEVENT;
0277 }
0278
0279 if (Verbosity())
0280 {
0281 m_central->identify();
0282 }
0283
0284 return Fun4AllReturnCodes::EVENT_OK;
0285 }
0286
0287 int CentralityReco::GetNodes(PHCompositeNode *topNode)
0288 {
0289 if (Verbosity() > 1)
0290 {
0291 std::cout << __FILE__ << " :: " << __FUNCTION__ << " :: " << __LINE__ << std::endl;
0292 }
0293
0294 m_mb_info = findNode::getClass<MinimumBiasInfo>(topNode, "MinimumBiasInfo");
0295
0296 if (!m_mb_info)
0297 {
0298 std::cout << "no mb_info node " << std::endl;
0299 return Fun4AllReturnCodes::ABORTRUN;
0300 }
0301
0302 m_global_vertex_map = findNode::getClass<GlobalVertexMap>(topNode, "GlobalVertexMap");
0303
0304 if (!m_global_vertex_map)
0305 {
0306 std::cout << "no vertex map node " << std::endl;
0307 return Fun4AllReturnCodes::EVENT_OK;
0308 }
0309
0310 m_central = findNode::getClass<CentralityInfo>(topNode, "CentralityInfo");
0311
0312 if (!m_central)
0313 {
0314 std::cout << "no centrality node " << std::endl;
0315 return Fun4AllReturnCodes::ABORTRUN;
0316 }
0317
0318 m_mbd_container = findNode::getClass<MbdPmtContainer>(topNode, "MbdPmtContainer");
0319
0320 if (!m_mbd_container)
0321 {
0322 std::cout << "no MBD out node " << std::endl;
0323 return Fun4AllReturnCodes::ABORTRUN;
0324 }
0325
0326
0327 m_mbd_out = findNode::getClass<MbdOut>(topNode, "MbdOut");
0328 if (Verbosity())
0329 {
0330 std::cout << "Getting MBD Out" << std::endl;
0331 }
0332
0333 if (!m_mbd_out)
0334 {
0335 std::cout << "no MBD out node " << std::endl;
0336 return Fun4AllReturnCodes::ABORTRUN;
0337 }
0338
0339 return Fun4AllReturnCodes::EVENT_OK;
0340 }
0341
0342 void CentralityReco::CreateNodes(PHCompositeNode *topNode)
0343 {
0344 if (Verbosity())
0345 {
0346 std::cout << __FILE__ << " :: " << __FUNCTION__ << " :: " << __LINE__ << std::endl;
0347 }
0348
0349 PHNodeIterator iter(topNode);
0350 PHCompositeNode *dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
0351 if (!dstNode)
0352 {
0353 std::cout << PHWHERE << "DST Node missing doing nothing" << std::endl;
0354 }
0355
0356 PHNodeIterator dstIter(dstNode);
0357
0358 PHCompositeNode *detNode = dynamic_cast<PHCompositeNode *>(dstIter.findFirst("PHCompositeNode", "GLOBAL"));
0359 if (!detNode)
0360 {
0361 std::cout << PHWHERE << "Detector Node missing, making one" << std::endl;
0362 detNode = new PHCompositeNode("GLOBAL");
0363 dstNode->addNode(detNode);
0364 }
0365
0366 CentralityInfo *central = new CentralityInfov2();
0367
0368 PHIODataNode<PHObject> *centralityNode = new PHIODataNode<PHObject>(central, "CentralityInfo", "PHObject");
0369 detNode->addNode(centralityNode);
0370
0371 return;
0372 }
0373
0374 float CentralityReco::getVertexScale()
0375 {
0376
0377
0378 float mbd_vertex = m_mbd_out->get_zvtx();
0379
0380 for (auto v_range_scale : m_vertex_scales)
0381 {
0382 auto v_range = v_range_scale.first;
0383 if (Verbosity())
0384 {
0385 std::cout << "vertexrange : "<<v_range.first<<"-"<<v_range.second << std::endl;
0386 }
0387
0388 if (mbd_vertex > v_range.first && mbd_vertex <= v_range.second)
0389 {
0390 return v_range_scale.second;
0391 }
0392 }
0393 return 0;
0394 }
0395
