Skip to content
Draft
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
249 changes: 166 additions & 83 deletions offline/packages/jetbackground/DetermineTowerBackground.cc
Original file line number Diff line number Diff line change
Expand Up @@ -14,6 +14,10 @@
#include <eventplaneinfo/Eventplaneinfo.h>
#include <eventplaneinfo/EventplaneinfoMap.h>

#include <centrality/CentralityInfo.h>
#include <ffamodules/CDBInterface.h>
#include <cdbobjects/CDBTTree.h>

#include <jetbase/Jet.h>
#include <jetbase/JetContainer.h>

Expand Down Expand Up @@ -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)
{

Expand Down Expand Up @@ -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<PHG4TruthInfoContainer>(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<EventplaneinfoMap>(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)
Expand Down Expand Up @@ -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<PHG4TruthInfoContainer>(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<EventplaneinfoMap>(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))
{
Expand Down Expand Up @@ -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<CentralityInfo>(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
Expand Down
14 changes: 13 additions & 1 deletion offline/packages/jetbackground/DetermineTowerBackground.h
Original file line number Diff line number Diff line change
Expand Up @@ -13,6 +13,7 @@
#include <jetbase/Jet.h>
#include <string>
#include <vector>
#include <array>

// forward declarations
class PHCompositeNode;
Expand All @@ -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; };
Expand All @@ -55,6 +60,13 @@ class DetermineTowerBackground : public SubsysReco
int CreateNode(PHCompositeNode *topNode);
void FillNode(PHCompositeNode *topNode);

int LoadCalibrations();

std::vector<float> _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};
Expand Down
3 changes: 3 additions & 0 deletions offline/packages/jetbackground/Makefile.am
Original file line number Diff line number Diff line change
Expand Up @@ -24,6 +24,8 @@ libjetbackground_la_LDFLAGS = \
libjetbackground_la_LIBADD = \
libjetbackground_io.la \
-lcalo_io \
-lcentrality_io \
-lcdbobjects \
-lConstituentSubtractor \
-leventplaneinfo_io \
-lglobalvertex \
Expand All @@ -33,6 +35,7 @@ libjetbackground_la_LIBADD = \
-lphg4hit \
-lphparameter \
-lqautils \
-lffamodules \
-lSubsysReco

pkginclude_HEADERS = \
Expand Down