src/TowerRho/DetermineTowerRho.cc

0001 #include "DetermineTowerRho.h"
0002
0003 #include "TowerRho.h"
0004 #include "TowerRhov1.h"
0005
0006 #include <jetbase/JetInput.h>
0007 #include <jetbase/Jet.h>
0008
0009 #include <fun4all/Fun4AllReturnCodes.h>
0010 #include <fun4all/SubsysReco.h>
0011
0012 #include <phool/PHCompositeNode.h>
0013 #include <phool/PHIODataNode.h>
0014 #include <phool/PHNode.h>
0015 #include <phool/PHNodeIterator.h>
0016 #include <phool/PHObject.h>
0017 #include <phool/getClass.h>
0018 #include <phool/phool.h>
0019
0020 // standard includes
0021 #include <cstdlib> // for exit
0022 #include <iostream>
0023 #include <vector>
0024 #include <string>
0025 #include <cmath>
0026 #include <algorithm>
0027
0028 // fastjet includes
0029 #include <fastjet/AreaDefinition.hh>
0030 #include <fastjet/ClusterSequenceArea.hh>
0031 #include <fastjet/Selector.hh>
0032 #include <fastjet/tools/BackgroundEstimatorBase.hh>
0033 #include <fastjet/tools/JetMedianBackgroundEstimator.hh>
0034 #include <fastjet/ClusterSequence.hh>
0035 #include <fastjet/JetDefinition.hh>
0036 #include <fastjet/PseudoJet.hh>
0037
0038 DetermineTowerRho::DetermineTowerRho(const std::string &name)
0039   : SubsysReco(name)
0040 {
0041 }
0042
0043 DetermineTowerRho::~DetermineTowerRho()
0044 {
0045   for (auto & _input : _inputs)
0046   {
0047     delete _input;
0048   }
0049   _inputs.clear();
0050   _outputs.clear();
0051 }
0052
0053 int DetermineTowerRho::InitRun(PHCompositeNode *topNode)
0054 {
0055   if(Verbosity() > 0) print_settings(std::cout);
0056   m_fj_algo_name = get_fj_algo_name();
0057   return CreateNode(topNode);
0058 }
0059
0060
0061 int DetermineTowerRho::process_event(PHCompositeNode *topNode)
0062 {
0063
0064     if (Verbosity() > 1) std::cout << "DetermineTowerRho::process_event -- entered" << std::endl;
0065
0066     if(m_abs_jet_eta_range == 5.0)
0067     {
0068       m_abs_jet_eta_range = m_abs_tower_eta_range - m_par;
0069
0070     }
0071     if(m_abs_ghost_eta == 5.0)
0072     {
0073       m_abs_ghost_eta = m_abs_tower_eta_range;
0074     }
0075
0076     // get fastjet input from tower nodes
0077     std::vector<fastjet::PseudoJet> calo_pseudojets = get_pseudojets(topNode);
0078
0079     // fastjet definitions
0080     fastjet::JetAlgorithm _fj_algo = get_fj_algo(); // get fastjet algorithm
0081
0082     // jet definition
0083     fastjet::JetDefinition jet_def(_fj_algo,
0084       m_par,
0085       fastjet::E_scheme,
0086       fastjet::Best);
0087
0088     // area definition
0089     fastjet::AreaDefinition area_def(fastjet::active_area_explicit_ghosts,
0090                   fastjet::GhostedAreaSpec(m_abs_ghost_eta, 1, m_ghost_area));
0091
0092
0093     fastjet::Selector jet_selector = (
0094       (!fastjet::SelectorNHardest(m_omit_nhardest)) *
0095       (fastjet::SelectorAbsEtaMax(m_abs_jet_eta_range) && fastjet::SelectorPtMin(m_jet_min_pT))
0096     );
0097
0098     // Not pure ghost function
0099     fastjet::Selector not_pure_ghost = (!fastjet::SelectorIsPureGhost());
0100
0101     for (unsigned int ipos = 0; ipos<_rho_methods.size(); ipos++)
0102     {
0103
0104       TowerRho::Method _rho_method = _rho_methods.at(ipos);
0105       float rho = 0;
0106       float sigma = 0;
0107
0108       if (_rho_method == TowerRho::Method::AREA)
0109       {
0110           fastjet::JetMedianBackgroundEstimator bge {jet_selector, jet_def, area_def};
0111           bge.set_particles(calo_pseudojets);
0112           rho = bge.rho();
0113           sigma = bge.sigma();
0114       }
0115       else if (_rho_method == TowerRho::Method::MULT)
0116       {
0117
0118               // reconstruct the background jets
0119             fastjet::ClusterSequenceArea cs(calo_pseudojets, jet_def, area_def);
0120             std::vector<fastjet::PseudoJet> jets = fastjet::sorted_by_pt( jet_selector( cs.inclusive_jets() ) );
0121
0122             std::vector<float> pt_over_nConstituents;
0123             int nfj_jets = 0;
0124             int n_empty_jets = 0;
0125             float total_constituents = 0;
0126             for (auto jet : jets)
0127             {
0128               float jet_pt = jet.pt();
0129
0130               std::vector<fastjet::PseudoJet> constituents = not_pure_ghost(jet.constituents());
0131
0132               int jet_size = constituents.size();
0133
0134               total_constituents += jet_size;
0135               if(jet_size == 0)
0136               {
0137                 n_empty_jets++;
0138                 continue; // only consider jets with constituents
0139               }
0140               pt_over_nConstituents.push_back(jet_pt/jet_size);
0141               nfj_jets++;
0142             }
0143
0144             int total_jets = nfj_jets + n_empty_jets;
0145             float mean_N = total_constituents/total_jets;
0146             if (mean_N < 0) mean_N = 0;
0147
0148             float med_tmp, std_tmp;
0149             _median_stddev(pt_over_nConstituents, n_empty_jets, med_tmp, std_tmp);
0150             sigma = std_tmp*sqrt(mean_N);
0151             rho = med_tmp;
0152
0153       }
0154       else{
0155         rho = 0;
0156         sigma = 0;
0157       }
0158
0159       FillNode(topNode, ipos, rho, sigma);
0160     }
0161
0162     return Fun4AllReturnCodes::EVENT_OK;
0163 }
0164
0165 std::vector<fastjet::PseudoJet> DetermineTowerRho::get_pseudojets(PHCompositeNode *topNode)
0166 {
0167     // get fastjet input from tower nodes
0168     std::vector<fastjet::PseudoJet> calo_pseudojets;
0169     for (auto & _input : _inputs)
0170     {
0171       std::vector<Jet *> parts = _input->get_input(topNode);
0172       for (unsigned int i = 0; i < parts.size(); ++i)
0173       {
0174         auto& part = parts[i];
0175         float this_e = part->get_e();
0176         if (this_e == 0.) continue;
0177         float this_px = parts[i]->get_px();
0178         float this_py = parts[i]->get_py();
0179         float this_pz = parts[i]->get_pz();
0180
0181
0182         if(m_do_tower_cut)
0183         {
0184           if(this_e < m_tower_threshold) continue;
0185           if (this_e < 0)
0186           {
0187             // make energy = +1 MeV for purposes of clustering
0188             float e_ratio = 0.001 / this_e;
0189             this_e  = this_e * e_ratio;
0190             this_px = this_px * e_ratio;
0191             this_py = this_py * e_ratio;
0192             this_pz = this_pz * e_ratio;
0193           }
0194         }
0195         fastjet::PseudoJet pseudojet(this_px, this_py, this_pz, this_e);
0196
0197         if(pseudojet.perp() < m_tower_min_pT) continue;
0198         if(fabs(pseudojet.eta()) > m_abs_tower_eta_range) continue;
0199
0200         calo_pseudojets.push_back(pseudojet);
0201       }
0202       for (auto &p : parts) delete p;
0203       parts.clear();
0204
0205     }
0206
0207     return calo_pseudojets;
0208 }
0209
0210
0211 int DetermineTowerRho::CreateNode(PHCompositeNode *topNode)
0212 {
0213   PHNodeIterator iter(topNode);
0214
0215   // Looking for the DST node
0216   PHCompositeNode *dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
0217   if (!dstNode)
0218   {
0219     std::cout << PHWHERE << "DST Node missing, doing nothing." << std::endl;
0220     return Fun4AllReturnCodes::ABORTRUN;
0221   }
0222
0223   // store the jet background stuff under a sub-node directory
0224   PHCompositeNode *bkgNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "JETBACKGROUND"));
0225   if (!bkgNode)
0226   {
0227     bkgNode = new PHCompositeNode("JETBACKGROUND");
0228     dstNode->addNode(bkgNode);
0229   }
0230
0231   // create the TowerRho nodes
0232   for (auto &_output : _outputs)
0233   {
0234     TowerRho *towerrho = findNode::getClass<TowerRho>(topNode, _output);
0235     if (!towerrho)
0236     {
0237       towerrho = new TowerRhov1();
0238       PHIODataNode<PHObject> *rhoDataNode = new PHIODataNode<PHObject>(towerrho, _output, "PHObject");
0239       bkgNode->addNode(rhoDataNode);
0240     }
0241   }
0242
0243
0244   return Fun4AllReturnCodes::EVENT_OK;
0245 }
0246
0247 void DetermineTowerRho::FillNode(PHCompositeNode *topNode, unsigned int ipos, const float rho, const float sigma)
0248 {
0249   std::string _node_name = _outputs.at(ipos);
0250   TowerRho::Method rho_method = _rho_methods.at(ipos);
0251   TowerRho *towerbackground = findNode::getClass<TowerRho>(topNode, _node_name);
0252   if (!towerbackground)
0253   {
0254     std::cout << " ERROR -- can't find TowerRho node after it should have been created" << std::endl;
0255     return;
0256   }
0257   else
0258   {
0259     if(Verbosity() > 0) std::cout << "DetermineTowerRho::FillNode - Filling node " << _node_name << " with rho = " << rho << " and sigma = " << sigma << std::endl;
0260     towerbackground->set_rho(rho);
0261     towerbackground->set_sigma(sigma);
0262     towerbackground->set_method(rho_method);
0263   }
0264
0265   return;
0266 }
0267
0268 float DetermineTowerRho::_percentile(const std::vector<float> & sorted_vec,
0269         const float percentile,
0270         const float nempty) const
0271 {
0272   assert(percentile >= 0. && percentile <= 1.);
0273
0274   int njets = sorted_vec.size();
0275   if(njets == 0) return 0;
0276
0277   float total_njets = njets + nempty;
0278   float perc_pos = (total_njets)*percentile - nempty - 0.5;
0279
0280   float result;
0281   if(perc_pos >= 0 && njets > 1){
0282     int pindex = int(perc_pos);
0283     // avoid out of range
0284     if(pindex + 1 > njets-1){
0285       pindex = njets-2;
0286       perc_pos = njets - 1;
0287     }
0288
0289     result = sorted_vec.at(pindex)*(pindex+1-perc_pos)
0290     + sorted_vec.at(pindex+1)*(perc_pos-pindex);
0291   }
0292   else if(perc_pos > -0.5 && njets >=1){
0293     result = sorted_vec.at(0);
0294   }
0295   else{
0296     result = 0;
0297   }
0298
0299   return result;
0300
0301 }
0302
0303 void DetermineTowerRho::_median_stddev(const std::vector<float> & vec,
0304               float n_empty_jets,
0305               float & median,
0306               float & std_dev) const
0307 {
0308   if(vec.size() == 0){
0309     median = 0;
0310     std_dev = 0;
0311     return;
0312   }
0313
0314   std::vector<float> sorted_vec = vec;
0315   std::sort(sorted_vec.begin(), sorted_vec.end());
0316
0317   int njets = sorted_vec.size();
0318   if(n_empty_jets > njets/4.0){
0319     std::cout << "WARNING: n_empty_jets = " << n_empty_jets << " is too large, setting to " << njets/4.0 << std::endl;
0320     n_empty_jets = njets/4.0;
0321   }
0322
0323
0324   float posn[2] = {0.5, (1.0-0.6827)/2.0};
0325   float res[2];
0326   for (int i = 0; i < 2; i++) {
0327     res[i] = _percentile(sorted_vec, posn[i], n_empty_jets);
0328   }
0329
0330   median = res[0];
0331   std_dev = res[0] - res[1];
0332
0333 }
0334
0335
0336 fastjet::JetAlgorithm DetermineTowerRho::get_fj_algo() {
0337
0338   if (_algo == Jet::ANTIKT)
0339   {
0340     return fastjet::antikt_algorithm;
0341   }
0342   else if (_algo == Jet::KT)
0343   {
0344     return fastjet::kt_algorithm;
0345   }
0346   else if (_algo == Jet::CAMBRIDGE)
0347   {
0348     return fastjet::cambridge_algorithm;
0349   }
0350   else
0351   {
0352     std::cout << PHWHERE << " jet algorithm not recognized, using default (kt)." << std::endl;
0353     return fastjet::kt_algorithm;
0354   }
0355
0356 }
0357
0358 std::string DetermineTowerRho::get_fj_algo_name(){
0359   if (_algo == Jet::ANTIKT){
0360     return "antikt_algorithm";
0361   }
0362   else if (_algo == Jet::KT){
0363     return "kt_algorithm";
0364   }
0365   else if (_algo == Jet::CAMBRIDGE){
0366     return "cambridge_algorithm";
0367   }
0368   else{
0369     std::cout << PHWHERE << " jet algorithm not recognized, using default (kt)." << std::endl;
0370     return "kt_algorithm";
0371   }
0372 }
0373
0374 void DetermineTowerRho::print_settings(std::ostream& os) const
0375 {
0376
0377   os << "DetermineTowerRho settings: " << std::endl;
0378   os << "Methods: ";
0379   for (auto & _rho_method : _rho_methods){
0380     os << TowerRhov1::get_method_string(_rho_method) << ", ";
0381   }
0382   os << std::endl;
0383   os << "Inputs:";
0384   for (auto & _input : _inputs) _input->identify();
0385   os << "Outputs: ";
0386   for (auto & _output : _outputs) os << _output << ", ";
0387   os << std::endl;
0388   os << "Using jet algo: " << m_fj_algo_name << " with R = " << m_par << std::endl;
0389   os << "Tower eta range: " << m_abs_tower_eta_range << std::endl;
0390   os << "Jet eta range: " << m_abs_jet_eta_range << std::endl;
0391   os << "Ghost area: " << m_ghost_area << std::endl;
0392   os << "Omit n hardest: " << m_omit_nhardest << std::endl;
0393   os << "Tower min pT: " << m_tower_min_pT << std::endl;
0394   os << "Jet min pT: " << m_jet_min_pT << std::endl;
0395   os << "Ghost eta: " << m_abs_ghost_eta << std::endl;
0396   os << "-----------------------------------" << std::endl;
0397   return;
0398
0399 }