packages/jetbackground/DetermineEventRho.cc

0001 #include "DetermineEventRho.h"
0002
0003 #include "EventRhov1.h"
0004
0005 #include <fun4all/SubsysReco.h>
0006 #include <jetbase/JetInput.h>
0007
0008 #include <fun4all/Fun4AllReturnCodes.h>
0009
0010 #include <phool/PHCompositeNode.h>
0011 #include <phool/PHIODataNode.h>
0012 #include <phool/PHNode.h>  // for PHNode
0013 #include <phool/PHNodeIterator.h>
0014 #include <phool/PHObject.h>
0015 #include <phool/getClass.h>
0016 #include <phool/phool.h>
0017
0018 // fastjet includes
0019 #include <fastjet/AreaDefinition.hh>
0020 #include <fastjet/ClusterSequence.hh>
0021 #include <fastjet/ClusterSequenceArea.hh>
0022 #include <fastjet/GhostedAreaSpec.hh>  // for GhostedAreaSpec
0023 #include <fastjet/JetDefinition.hh>
0024 #include <fastjet/PseudoJet.hh>
0025 #include <fastjet/Selector.hh>
0026
0027 // standard includes
0028 #include <algorithm>
0029 #include <cmath>
0030 #include <cstdlib>
0031 #include <iostream>
0032 #include <sstream>  // for basic_ostringstream
0033 #include <string>
0034 #include <vector>
0035
0036 DetermineEventRho::DetermineEventRho(const std::string &name)
0037   : SubsysReco(name)
0038 {
0039   // silence output from fastjet
0040   fastjet::ClusterSequence const clusseq;
0041   if (Verbosity() > 0)
0042   {
0043     fastjet::ClusterSequence::print_banner();
0044   }
0045   else
0046   {
0047     std::ostringstream nullstream;
0048     fastjet::ClusterSequence::set_fastjet_banner_stream(&nullstream);
0049     fastjet::ClusterSequence::print_banner();
0050     fastjet::ClusterSequence::set_fastjet_banner_stream(&std::cout);
0051   }
0052 }
0053
0054 DetermineEventRho::~DetermineEventRho()
0055 {
0056   // clean up memory
0057   for (auto &input : m_inputs)
0058   {
0059     delete input;
0060   }
0061   m_inputs.clear();
0062   m_output_nodes.clear();
0063   m_rho_methods.clear();
0064 }
0065
0066 int DetermineEventRho::InitRun(PHCompositeNode *topNode)
0067 {
0068   // set jet_eta_range
0069   if (m_abs_jet_eta_range < 0)
0070   {
0071     m_abs_jet_eta_range = m_abs_input_eta_range - m_par;
0072   }
0073   if (Verbosity() > 0)
0074   {
0075     print_settings();
0076   }
0077
0078   return CreateNodes(topNode);
0079 }
0080
0081 int DetermineEventRho::process_event(PHCompositeNode *topNode)
0082 {
0083   if (Verbosity() > 1)
0084   {
0085     std::cout << "DetermineEventRho::process_event -- entered" << std::endl;
0086   }
0087
0088   std::vector<Jet *> particles{};
0089   for (auto &input : m_inputs)
0090   {
0091     std::vector<Jet *> const parts = input->get_input(topNode);
0092     for (const auto &part : parts)
0093     {
0094       particles.push_back(part);
0095       particles.back()->set_id(particles.size() - 1);  // unique ids ensured
0096     }
0097   }
0098
0099   std::vector<fastjet::PseudoJet> calo_pseudojets{};
0100   for (unsigned int ipart = 0; ipart < particles.size(); ++ipart)
0101   {
0102     float this_e = particles[ipart]->get_e();
0103     if (this_e == 0.)
0104     {
0105       continue;
0106     }  // skip zero energy particles
0107     float this_px = particles[ipart]->get_px();
0108     float this_py = particles[ipart]->get_py();
0109     float this_pz = particles[ipart]->get_pz();
0110
0111     if (this_e < 0)
0112     {  // make energy = +1 MeV for purposes of clustering
0113       float const e_ratio = 0.001 / this_e;
0114       this_e = this_e * e_ratio;
0115       this_px = this_px * e_ratio;
0116       this_py = this_py * e_ratio;
0117       this_pz = this_pz * e_ratio;
0118     }
0119
0120     fastjet::PseudoJet pseudojet(this_px, this_py, this_pz, this_e);
0121     // eta cut
0122     if (std::abs(pseudojet.eta()) > m_abs_input_eta_range)
0123     {
0124       continue;
0125     }
0126
0127     pseudojet.set_user_index(ipart);
0128     calo_pseudojets.push_back(pseudojet);
0129
0130   }  // end of loop over particles
0131
0132   // initialize the jet selector
0133   auto jet_selector = get_jet_selector();
0134   fastjet::JetDefinition *m_jet_def = nullptr;
0135   if (m_bkgd_jet_algo == Jet::ALGO::ANTIKT)
0136   {
0137     m_jet_def = new fastjet::JetDefinition(fastjet::antikt_algorithm, m_par, fastjet::E_scheme, fastjet::Best);
0138   }
0139   else if (m_bkgd_jet_algo == Jet::ALGO::KT)
0140   {
0141     m_jet_def = new fastjet::JetDefinition(fastjet::kt_algorithm, m_par, fastjet::E_scheme, fastjet::Best);
0142   }
0143   else if (m_bkgd_jet_algo == Jet::ALGO::CAMBRIDGE)
0144   {
0145     m_jet_def = new fastjet::JetDefinition(fastjet::cambridge_algorithm, m_par, fastjet::E_scheme, fastjet::Best);
0146   }
0147   else
0148   {
0149     std::cout << PHWHERE << " jet algorithm not recognized, using default (kt)." << std::endl;
0150     m_jet_def = new fastjet::JetDefinition(fastjet::kt_algorithm, m_par, fastjet::E_scheme, fastjet::Best);
0151   }
0152   for (unsigned int ipos = 0; ipos < m_rho_methods.size(); ipos++)
0153   {
0154     float rho = 0;
0155     float sigma = 0;
0156     auto rho_method = m_rho_methods.at(ipos);
0157
0158     auto *m_eventbackground = findNode::getClass<EventRho>(topNode, m_output_nodes.at(ipos));
0159     if (!m_eventbackground)
0160     {
0161       std::cout << PHWHERE << " EventRho node " << m_output_nodes.at(ipos) << " not found" << std::endl;
0162       continue;
0163     }
0164
0165     if (!m_jet_def)
0166     {
0167       std::cerr << PHWHERE << " jet definition not set" << std::endl;
0168       exit(1);
0169     }
0170
0171     if (rho_method == EventRho::Method::AREA)
0172     {
0173       fastjet::AreaDefinition const area_def(fastjet::active_area_explicit_ghosts,
0174                                              fastjet::GhostedAreaSpec(m_abs_input_eta_range, 1, m_ghost_area));
0175       auto *m_cluseq = new fastjet::ClusterSequenceArea(calo_pseudojets, *m_jet_def, area_def);
0176       auto fastjets = jet_selector(m_cluseq->inclusive_jets());
0177
0178       std::vector<float> pT_over_X{};
0179       float total_X = 0;
0180       float njets_used = 0.0;
0181       float empty_X = 0;
0182       float const njets_total = static_cast<float>(fastjets.size());
0183
0184       for (auto &fastjet : fastjets)
0185       {
0186         float const this_X = fastjet.area();
0187         // if (this_X <= 0 || this_X != this_X || fastjet.is_pure_ghost())
0188         if (this_X <= 0 || this_X != this_X)
0189         {
0190           if (Verbosity() > 2)
0191           {
0192             std::cout << PHWHERE << " ::WARNING: Discarding jet with zero area. Zero-area jets may be due to (i) too large a ghost area (ii) a jet being outside the ghost range (iii) the computation not being done using an appropriate algorithm (kt;C/A)." << std::endl;
0193           }
0194           if (!std::isnan(this_X))
0195           {
0196             empty_X += this_X;
0197           }
0198           continue;  // skip this jet
0199         }  // end of check on X
0200
0201         // add back the truth comp pT
0202         float px_sum = 0;
0203         float py_sum = 0;
0204         for (auto &comp : fastjet.constituents())
0205         {
0206           if (comp.is_pure_ghost())
0207           {  // skip pure ghosts
0208             continue;
0209           }
0210
0211           auto *p = particles[comp.user_index()];  // get original particle
0212           px_sum += p->get_px();
0213           py_sum += p->get_py();
0214
0215         }  // end of loop over fastjet comps
0216
0217         // calculate pT and pT/X
0218         float const this_pT_over_X = std::sqrt((px_sum * px_sum) + (py_sum * py_sum)) / this_X;
0219         pT_over_X.push_back(this_pT_over_X);
0220         total_X += this_X;
0221         njets_used += 1.0;
0222       }  // end of loop over fastjets
0223
0224       if (empty_X != 0.0)
0225       {
0226         if (Verbosity() > 0)
0227         {
0228           std::cerr << PHWHERE << " ::WARNING: Found " << empty_X << " empty jets with zero area. This may be due to (i) too large a ghost area (ii) a jet being outside the ghost range (iii) the computation not being done using an appropriate algorithm (kt;C/A)." << std::endl;
0229         }
0230         total_X += empty_X;
0231       }
0232
0233       float const n_empty_jets = njets_total - njets_used;
0234       float mean_X = (1.0 * total_X) / (njets_total);
0235       if (mean_X < 0)
0236       {
0237         std::cerr << PHWHERE << " mean_N < 0 , setting to 0" << std::endl;
0238         mean_X = 0;
0239       }
0240
0241       float tmp_med;
0242       float tmp_std;
0243       CalcMedianStd(pT_over_X, n_empty_jets, tmp_med, tmp_std);
0244
0245       sigma = std::sqrt(mean_X) * tmp_std;
0246       rho = tmp_med;
0247
0248       // clean up
0249       fastjets.clear();
0250       pT_over_X.clear();
0251       delete m_cluseq;
0252     }
0253     else if (rho_method == EventRho::Method::MULT)
0254     {
0255       // reconstruct the background jets
0256       auto *m_cluseq = new fastjet::ClusterSequence(calo_pseudojets, *m_jet_def);
0257       auto fastjets = jet_selector(m_cluseq->inclusive_jets());
0258
0259       std::vector<float> pt_over_nconst{};
0260       int total_constituents = 0;
0261
0262       for (auto &fastjet : fastjets)
0263       {
0264         auto comps = fastjet.constituents();
0265         if (!comps.empty())
0266         {
0267           float total_px = 0;
0268           float total_py = 0;
0269           for (auto &comp : comps)
0270           {
0271             auto *particle = particles[comp.user_index()];
0272             total_px += particle->get_px();
0273             total_py += particle->get_py();
0274             total_constituents++;
0275           }  // end of loop over constituents
0276           float const jet_avg_pt = (std::sqrt((total_px * total_px) + (total_py * total_py)) / (1.0 * comps.size()));
0277           pt_over_nconst.push_back(jet_avg_pt);
0278
0279         }  // end of if comps.size() > 0
0280       }  // end of loop over fastjets
0281
0282       float const n_empty_jets = 1.0 * (fastjets.size() - pt_over_nconst.size());
0283       float mean_N = (1.0 * total_constituents) / (1.0 * fastjets.size());
0284       if (mean_N < 0)
0285       {
0286         std::cerr << PHWHERE << " mean_N < 0 , setting to 0" << std::endl;
0287         mean_N = 0;
0288       }
0289
0290       float tmp_med;
0291       float tmp_std;
0292       CalcMedianStd(pt_over_nconst, 1.0 * n_empty_jets, tmp_med, tmp_std);
0293
0294       sigma = std::sqrt(mean_N) * tmp_std;
0295       rho = tmp_med;
0296
0297       // clean up
0298       fastjets.clear();
0299       pt_over_nconst.clear();
0300       delete m_cluseq;
0301     }
0302     else
0303     {
0304       std::cout << PHWHERE << " rho method not recognized" << std::endl;
0305     }
0306
0307     if (Verbosity() > 1)
0308     {
0309       std::cout << "DetermineEventRho::process_event - Filling node "
0310                 << m_output_nodes.at(ipos) << " with rho = " << rho
0311                 << " and sigma = " << sigma << std::endl;
0312     }
0313
0314     m_eventbackground->set_rho(rho);
0315     m_eventbackground->set_sigma(sigma);
0316     m_eventbackground->set_method(rho_method);
0317
0318   }  // end of loop over rho methods
0319
0320   delete m_jet_def;
0321   // clean up input vectors
0322   for (auto &part : particles)
0323   {
0324     delete part;
0325   }
0326   particles.clear();
0327   calo_pseudojets.clear();
0328
0329   return Fun4AllReturnCodes::EVENT_OK;
0330 }
0331
0332 void DetermineEventRho::add_method(EventRho::Method rho_method, std::string output_node)
0333 {
0334   // get method name ( also checks if method is valid )
0335   std::string const method_name = EventRhov1::get_method_string(rho_method);
0336
0337   // check if method already exists
0338   if (std::find(m_rho_methods.begin(), m_rho_methods.end(), rho_method) != m_rho_methods.end())
0339   {
0340     std::cout << PHWHERE << " method " << method_name << " already exists, skipping" << std::endl;
0341     return;
0342   }
0343
0344   m_rho_methods.push_back(rho_method);
0345
0346   // if no output name is specified, use default
0347   if (output_node.empty())
0348   {
0349     output_node = "EventRho_" + method_name;
0350   }
0351   m_output_nodes.push_back(output_node);
0352   return;
0353 }
0354
0355 int DetermineEventRho::CreateNodes(PHCompositeNode *topNode)
0356 {
0357   PHNodeIterator iter(topNode);  // Looking for the DST node
0358   auto *dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
0359   if (!dstNode)
0360   {
0361     std::cout << PHWHERE << "DST Node missing, doing nothing." << std::endl;
0362     return Fun4AllReturnCodes::ABORTRUN;
0363   }
0364   auto *bkgNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "JETBACKGROUND"));
0365   if (!bkgNode)
0366   {  // create the node if it does not exist
0367     bkgNode = new PHCompositeNode("JETBACKGROUND");
0368     dstNode->addNode(bkgNode);
0369   }
0370
0371   // create the EventRho nodes
0372   for (auto &output : m_output_nodes)
0373   {
0374     auto *eventrho = findNode::getClass<EventRho>(topNode, output);
0375     if (!eventrho)
0376     {
0377       eventrho = new EventRhov1();
0378       auto *rhoDataNode = new PHIODataNode<PHObject>(eventrho, output, "PHObject");
0379       bkgNode->addNode(rhoDataNode);
0380     }  // end of if eventrho
0381   }  // end of loop over output nodes
0382
0383   return Fun4AllReturnCodes::EVENT_OK;
0384 }
0385
0386 float DetermineEventRho::CalcPercentile(const std::vector<float> &sorted_vec, const float percentile, const float nempty)
0387 {
0388   float result = 0;
0389   if (!sorted_vec.empty())
0390   {
0391     int const njets = sorted_vec.size();
0392     float const total_njets = njets + nempty;
0393     float perc_pos = ((total_njets) *percentile) - nempty - 0.5;
0394
0395     if (perc_pos >= 0 && njets > 1)
0396     {
0397       int pindex = int(perc_pos);
0398       if (pindex + 1 > njets - 1)
0399       {  // avoid out of range
0400         pindex = njets - 2;
0401         perc_pos = njets - 1;
0402       }
0403       result = sorted_vec.at(pindex) * (pindex + 1 - perc_pos) + sorted_vec.at(pindex + 1) * (perc_pos - pindex);
0404     }
0405     else if (perc_pos > -0.5 && njets >= 1)
0406     {
0407       result = sorted_vec.at(0);
0408     }
0409     else
0410     {
0411       result = 0;
0412     }  // end of if criteria
0413
0414   }  // end of if sorted_vec.size() > 0
0415
0416   return result;
0417 }
0418
0419 void DetermineEventRho::CalcMedianStd(const std::vector<float> &vec, float n_empty_jets, float &median, float &std_dev)
0420 {
0421   median = 0;
0422   std_dev = 0;
0423   if (!vec.empty())
0424   {
0425     // sort the vector
0426     std::vector<float> sorted_vec = vec;
0427     std::sort(sorted_vec.begin(), sorted_vec.end());
0428
0429     int const njets = sorted_vec.size();
0430     if (n_empty_jets > njets / 4.0)
0431     {
0432       std::cout << "WARNING: n_empty_jets = " << n_empty_jets << " is too large, setting to " << njets / 4.0 << std::endl;
0433       n_empty_jets = njets / 4.0;
0434     }
0435
0436     float const posn[2] = {0.5, (1.0 - 0.6827) / 2.0};
0437     float res[2] = {0, 0};
0438     for (int i = 0; i < 2; i++)
0439     {
0440       res[i] = CalcPercentile(sorted_vec, posn[i], n_empty_jets);
0441     }
0442
0443     median = res[0];
0444     std_dev = res[0] - res[1];
0445     sorted_vec.clear();
0446   }  // end of if vec.size() > 0
0447
0448   return;
0449 }
0450
0451 fastjet::Selector DetermineEventRho::get_jet_selector() const
0452 {
0453   if (m_jet_min_pT != VOID_CUT)
0454   {
0455     return (!fastjet::SelectorNHardest(m_omit_nhardest)) * (fastjet::SelectorAbsEtaMax(m_abs_jet_eta_range)) && (fastjet::SelectorPtMin(m_jet_min_pT));
0456   }
0457   // default is no min jet pT
0458   return (!fastjet::SelectorNHardest(m_omit_nhardest)) * (fastjet::SelectorAbsEtaMax(m_abs_jet_eta_range));
0459 }
0460
0461 void DetermineEventRho::print_settings(std::ostream &os)
0462 {
0463   os << PHWHERE << "-----------------------------------" << std::endl;
0464   os << "Methods: ";
0465   for (auto rho_method : m_rho_methods)
0466   {
0467     os << EventRhov1::get_method_string(rho_method) << ", ";
0468   }
0469   os << std::endl;
0470
0471   os << "Inputs:";
0472   for (auto &input : m_inputs)
0473   {
0474     input->identify();
0475   }
0476
0477   os << "Outputs: ";
0478   for (const auto &output : m_output_nodes)
0479   {
0480     os << output << ", ";
0481   }
0482   os << std::endl;
0483
0484   os << "Using jet algo: ";
0485   if (m_bkgd_jet_algo == Jet::ALGO::ANTIKT)
0486   {
0487     os << "ANTIKT r=" << m_par;
0488   }
0489   else if (m_bkgd_jet_algo == Jet::ALGO::KT)
0490   {
0491     os << "KT r=" << m_par;
0492   }
0493   else if (m_bkgd_jet_algo == Jet::ALGO::CAMBRIDGE)
0494   {
0495     os << "CAMBRIDGE r=" << m_par;
0496   }
0497   os << std::endl;
0498
0499   os << "Tower eta range: " << m_abs_input_eta_range << std::endl;
0500   os << "Jet eta range: " << m_abs_jet_eta_range << std::endl;
0501   os << "Ghost area: " << m_ghost_area << std::endl;
0502   os << "Ghost eta: " << m_abs_input_eta_range << std::endl;
0503   os << "Omit n hardest: " << m_omit_nhardest << std::endl;
0504   if (m_jet_min_pT != VOID_CUT)
0505   {
0506     os << "Jet min pT: " << m_jet_min_pT << std::endl;
0507   }
0508   os << "-----------------------------------" << std::endl;
0509   return;
0510 }