diff --git a/offline/packages/jetbackground/DetermineTowerBackground.cc b/offline/packages/jetbackground/DetermineTowerBackground.cc index 9823d337f6..7f021f4f39 100644 --- a/offline/packages/jetbackground/DetermineTowerBackground.cc +++ b/offline/packages/jetbackground/DetermineTowerBackground.cc @@ -14,6 +14,10 @@ #include #include +#include +#include +#include + #include #include @@ -52,9 +56,62 @@ DetermineTowerBackground::DetermineTowerBackground(const std::string &name) int DetermineTowerBackground::InitRun(PHCompositeNode *topNode) { + if (_do_flow == 4) + { + if (Verbosity()) + { + std::cout << "Loading the average calo v2" << std::endl; + } + if (!LoadCalibrations()) + { + std::cout << "Load calibrations failed." << std::endl; + return Fun4AllReturnCodes::ABORTRUN; + } + + } + return CreateNode(topNode); } +int DetermineTowerBackground::LoadCalibrations() +{ + + CDBTTree *cdbtree_calo_v2 = nullptr; + + std::string calibdir; + if (m_overwrite_average_calo_v2) + { + calibdir = m_overwrite_average_calo_v2_path; + } + else + { + calibdir = CDBInterface::instance()->getUrl(m_calibName); + } + + if (calibdir.empty()) + { + std::cout << "Could not find and load histograms for EMCAL LUTs! defaulting to the identity table!" << std::endl; + exit(-1); + } + else + { + cdbtree_calo_v2 = new CDBTTree(calibdir); + } + + cdbtree_calo_v2->LoadCalibrations(); + + _CENTRALITY_V2.assign(100,0); + + for (int icent = 0; icent < 100; icent++) + { + _CENTRALITY_V2[icent] = cdbtree_calo_v2->GetFloatValue(icent, "jet_calo_v2"); + } + + delete cdbtree_calo_v2; + + return Fun4AllReturnCodes::EVENT_OK; +} + int DetermineTowerBackground::process_event(PHCompositeNode *topNode) { @@ -481,7 +538,92 @@ int DetermineTowerBackground::process_event(PHCompositeNode *topNode) } } - if ( _do_flow >= 1 ) + + // Get psi + if (_do_flow == 2) + { // HIJING truth flow extraction + PHG4TruthInfoContainer *truthinfo = findNode::getClass(topNode, "G4TruthInfo"); + + if (!truthinfo) + { + std::cout << "DetermineTowerBackground::process_event: FATAL , G4TruthInfo does not exist , cannot extract truth flow with do_flow = " << _do_flow << std::endl; + return -1; + } + + PHG4TruthInfoContainer::Range range = truthinfo->GetPrimaryParticleRange(); + + float Hijing_Qx = 0; + float Hijing_Qy = 0; + + for (PHG4TruthInfoContainer::ConstIterator iter = range.first; iter != range.second; ++iter) + { + PHG4Particle *g4particle = iter->second; + + if (truthinfo->isEmbeded(g4particle->get_track_id()) != 0) + { + continue; + } + + TLorentzVector t; + t.SetPxPyPzE(g4particle->get_px(), g4particle->get_py(), g4particle->get_pz(), g4particle->get_e()); + + float truth_pt = t.Pt(); + if (truth_pt < 0.4) + { + continue; + } + float truth_eta = t.Eta(); + if (std::fabs(truth_eta) > 1.1) + { + continue; + } + float truth_phi = t.Phi(); + int truth_pid = g4particle->get_pid(); + + if (Verbosity() > 10) + { + std::cout << "DetermineTowerBackground::process_event: determining truth flow, using particle w/ pt / eta / phi " << truth_pt << " / " << truth_eta << " / " << truth_phi << " , embed / PID = " << truthinfo->isEmbeded(g4particle->get_track_id()) << " / " << truth_pid << std::endl; + } + + Hijing_Qx += truth_pt * std::cos(2 * truth_phi); + Hijing_Qy += truth_pt * std::sin(2 * truth_phi); + } + + _Psi2 = std::atan2(Hijing_Qy, Hijing_Qx) / 2.0; + + if (Verbosity() > 0) + { + std::cout << "DetermineTowerBackground::process_event: flow extracted from Hijing truth particles, setting Psi2 = " << _Psi2 << " ( " << _Psi2 / M_PI << " * pi ) " << std::endl; + } + } + else if (_do_flow == 3 || _do_flow == 4) + { // sEPD event plane extraction + // get event plane map + EventplaneinfoMap *epmap = findNode::getClass(topNode, "EventplaneinfoMap"); + if (!epmap) + { + std::cout << "DetermineTowerBackground::process_event: FATAL, EventplaneinfoMap does not exist, cannot extract sEPD flow with do_flow = " << _do_flow << std::endl; + exit(-1); + } + if (!(epmap->empty())) + { + auto *EPDNS = epmap->get(EventplaneinfoMap::sEPDNS); + _Psi2 = EPDNS->get_shifted_psi(2); + } + else + { + _is_flow_failure = true; + _Psi2 = 0; + } + + if (Verbosity() > 0) + { + std::cout << "DetermineTowerBackground::process_event: flow extracted from sEPD, setting Psi2 = " << _Psi2 << " ( " << _Psi2 / M_PI << " * pi ) " << std::endl; + } + + } + + if ( _do_flow >= 1 && _do_flow < 4) { if (Verbosity() > 0) @@ -754,88 +896,6 @@ int DetermineTowerBackground::process_event(PHCompositeNode *topNode) { // Calo event plane _Psi2 = std::atan2(Q_y, Q_x) / 2.0; } - else if (_do_flow == 2) - { // HIJING truth flow extraction - PHG4TruthInfoContainer *truthinfo = findNode::getClass(topNode, "G4TruthInfo"); - - if (!truthinfo) - { - std::cout << "DetermineTowerBackground::process_event: FATAL , G4TruthInfo does not exist , cannot extract truth flow with do_flow = " << _do_flow << std::endl; - return -1; - } - - PHG4TruthInfoContainer::Range range = truthinfo->GetPrimaryParticleRange(); - - float Hijing_Qx = 0; - float Hijing_Qy = 0; - - for (PHG4TruthInfoContainer::ConstIterator iter = range.first; iter != range.second; ++iter) - { - PHG4Particle *g4particle = iter->second; - - if (truthinfo->isEmbeded(g4particle->get_track_id()) != 0) - { - continue; - } - - TLorentzVector t; - t.SetPxPyPzE(g4particle->get_px(), g4particle->get_py(), g4particle->get_pz(), g4particle->get_e()); - - float truth_pt = t.Pt(); - if (truth_pt < 0.4) - { - continue; - } - float truth_eta = t.Eta(); - if (std::fabs(truth_eta) > 1.1) - { - continue; - } - float truth_phi = t.Phi(); - int truth_pid = g4particle->get_pid(); - - if (Verbosity() > 10) - { - std::cout << "DetermineTowerBackground::process_event: determining truth flow, using particle w/ pt / eta / phi " << truth_pt << " / " << truth_eta << " / " << truth_phi << " , embed / PID = " << truthinfo->isEmbeded(g4particle->get_track_id()) << " / " << truth_pid << std::endl; - } - - Hijing_Qx += truth_pt * std::cos(2 * truth_phi); - Hijing_Qy += truth_pt * std::sin(2 * truth_phi); - } - - _Psi2 = std::atan2(Hijing_Qy, Hijing_Qx) / 2.0; - - if (Verbosity() > 0) - { - std::cout << "DetermineTowerBackground::process_event: flow extracted from Hijing truth particles, setting Psi2 = " << _Psi2 << " ( " << _Psi2 / M_PI << " * pi ) " << std::endl; - } - } - else if (_do_flow == 3) - { // sEPD event plane extraction - // get event plane map - EventplaneinfoMap *epmap = findNode::getClass(topNode, "EventplaneinfoMap"); - if (!epmap) - { - std::cout << "DetermineTowerBackground::process_event: FATAL, EventplaneinfoMap does not exist, cannot extract sEPD flow with do_flow = " << _do_flow << std::endl; - exit(-1); - } - if (!(epmap->empty())) - { - auto *EPDNS = epmap->get(EventplaneinfoMap::sEPDNS); - _Psi2 = EPDNS->get_shifted_psi(2); - } - else - { - _is_flow_failure = true; - _Psi2 = 0; - } - - if (Verbosity() > 0) - { - std::cout << "DetermineTowerBackground::process_event: flow extracted from sEPD, setting Psi2 = " << _Psi2 << " ( " << _Psi2 / M_PI << " * pi ) " << std::endl; - } - - } if (std::isnan(_Psi2) || std::isinf(_Psi2)) { @@ -890,7 +950,30 @@ int DetermineTowerBackground::process_event(PHCompositeNode *topNode) std::cout << "DetermineTowerBackground::process_event: flow extraction successful, Psi2 = " << _Psi2 << " ( " << _Psi2 / M_PI << " * pi ) , v2 = " << _v2 << std::endl; } } // if do flow + else if (_do_flow == 4) + { + CentralityInfo *centinfo = findNode::getClass(topNode, "CentralityInfo"); + + if (!centinfo) + { + std::cout << "DetermineTowerBackground::process_event: FATAL, CentralityInfo does not exist, cannot extract centrality with do_flow = " << _do_flow << std::endl; + exit(-1); + } + + int centrality_bin = centinfo->get_centrality_bin(CentralityInfo::PROP::mbd_NS); + + if (centrality_bin > 0 && centrality_bin < 95) + { + _v2 = _CENTRALITY_V2[centrality_bin]; + } + else + { + _v2 = 0; + _is_flow_failure = true; + _Psi2 = 0; + } + } // now calculate energy densities... _nTowers = 0; // store how many towers were used to determine bkg diff --git a/offline/packages/jetbackground/DetermineTowerBackground.h b/offline/packages/jetbackground/DetermineTowerBackground.h index a8a6d0209c..ed9e34bf9a 100644 --- a/offline/packages/jetbackground/DetermineTowerBackground.h +++ b/offline/packages/jetbackground/DetermineTowerBackground.h @@ -13,6 +13,7 @@ #include #include #include +#include // forward declarations class PHCompositeNode; @@ -37,7 +38,11 @@ class DetermineTowerBackground : public SubsysReco void SetBackgroundOutputName(const std::string &name) { _backgroundName = name; } void SetSeedType(int seed_type) { _seed_type = seed_type; } void SetFlow(int do_flow) { _do_flow = do_flow; }; - + void SetOverwriteCaloV2(std::string &url) + { + m_overwrite_average_calo_v2 = true; + m_overwrite_average_calo_v2_path = url; + } void SetSeedJetD(float D) { _seed_jet_D = D; }; void SetSeedJetPt(float pt) { _seed_jet_pt = pt; }; void SetSeedMaxConst(float max_const) { _seed_max_const = max_const; }; @@ -55,6 +60,13 @@ class DetermineTowerBackground : public SubsysReco int CreateNode(PHCompositeNode *topNode); void FillNode(PHCompositeNode *topNode); + int LoadCalibrations(); + + std::vector _CENTRALITY_V2; + std::string m_calibName = "JET_AVERAGE_CALO_V2_SEPD_PSI2"; + bool m_overwrite_average_calo_v2{false}; + std::string m_overwrite_average_calo_v2_path; + int _do_flow{0}; float _v2{0}; float _Psi2{0}; diff --git a/offline/packages/jetbackground/Makefile.am b/offline/packages/jetbackground/Makefile.am index 6fe9740a4e..ad45676fbe 100644 --- a/offline/packages/jetbackground/Makefile.am +++ b/offline/packages/jetbackground/Makefile.am @@ -24,6 +24,8 @@ libjetbackground_la_LDFLAGS = \ libjetbackground_la_LIBADD = \ libjetbackground_io.la \ -lcalo_io \ + -lcentrality_io \ + -lcdbobjects \ -lConstituentSubtractor \ -leventplaneinfo_io \ -lglobalvertex \ @@ -33,6 +35,7 @@ libjetbackground_la_LIBADD = \ -lphg4hit \ -lphparameter \ -lqautils \ + -lffamodules \ -lSubsysReco pkginclude_HEADERS = \