packages/jetbackground/DetermineTowerRho.cc

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