Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-08-06 08:17:45 

0001 #include "FastJetAlgo.h"
0002 
0003 #include "Jet.h"
0004 #include "JetContainer.h"
0005 #include "Jetv2.h"
0006 
0007 #include <phool/phool.h>
0008 
0009 #include <TClonesArray.h>
0010 #include <TSystem.h>
0011 
0012 // fastjet includes
0013 #include <fastjet/AreaDefinition.hh>
0014 #include <fastjet/ClusterSequence.hh>
0015 #include <fastjet/ClusterSequenceArea.hh>
0016 #include <fastjet/FunctionOfPseudoJet.hh>  // for FunctionOfPse...
0017 #include <fastjet/JetDefinition.hh>
0018 #include <fastjet/PseudoJet.hh>
0019 #include <fastjet/Selector.hh>
0020 #include <fastjet/contrib/ConstituentSubtractor.hh>
0021 #include <fastjet/contrib/RecursiveSymmetryCutBase.hh>  // for RecursiveSymm...
0022 #include <fastjet/contrib/SoftDrop.hh>
0023 #include <fastjet/tools/BackgroundEstimatorBase.hh>
0024 #include <fastjet/tools/GridMedianBackgroundEstimator.hh>
0025 #include <fastjet/tools/JetMedianBackgroundEstimator.hh>
0026 
0027 #include <boost/format.hpp>
0028 
0029 // standard includes
0030 #include <cassert>
0031 #include <cmath>  // for isfinite
0032 #include <fstream>
0033 #include <iostream>
0034 #include <map>      // for _Rb_tree_iterator
0035 #include <memory>   // for allocator_traits<>::value_type
0036 #include <string>   // for operator<<
0037 #include <utility>  // for pair
0038 #include <vector>
0039 
0040 FastJetAlgo::FastJetAlgo(const FastJetOptions& options)
0041   : m_opt{options}
0042 {
0043   fastjet::ClusterSequence clusseq;
0044   if (m_opt.verbosity > 0)
0045   {
0046     fastjet::ClusterSequence::print_banner();
0047   }
0048   else
0049   {
0050     std::ostringstream nullstream;
0051     fastjet::ClusterSequence::set_fastjet_banner_stream(&nullstream);
0052     fastjet::ClusterSequence::print_banner();
0053     fastjet::ClusterSequence::set_fastjet_banner_stream(&std::cout);
0054   }
0055 }
0056 
0057 void FastJetAlgo::identify(std::ostream& os)
0058 {
0059   os << "   FastJetAlgo: ";
0060   if (m_opt.algo == Jet::ANTIKT)
0061   {
0062     os << "ANTIKT r=" << m_opt.jet_R;
0063   }
0064   else if (m_opt.algo == Jet::KT)
0065   {
0066     os << "KT r=" << m_opt.jet_R;
0067   }
0068   else if (m_opt.algo == Jet::CAMBRIDGE)
0069   {
0070     os << "CAMBRIDGE r=" << m_opt.jet_R;
0071   }
0072   os << std::endl;
0073 
0074   m_opt.print(os);
0075 }
0076 
0077 fastjet::JetDefinition FastJetAlgo::get_fastjet_definition() const
0078 {
0079   if (m_opt.algo == Jet::ANTIKT)
0080   {
0081     return fastjet::JetDefinition(fastjet::antikt_algorithm, m_opt.jet_R, fastjet::E_scheme, fastjet::Best);
0082   }
0083   if (m_opt.algo == Jet::KT)
0084   {
0085     return fastjet::JetDefinition(fastjet::kt_algorithm, m_opt.jet_R, fastjet::E_scheme, fastjet::Best);
0086   }
0087   if (m_opt.algo == Jet::CAMBRIDGE)
0088   {
0089     return fastjet::JetDefinition(fastjet::cambridge_algorithm, m_opt.jet_R, fastjet::E_scheme, fastjet::Best);
0090   }
0091 
0092   std::cout << PHWHERE << std::endl;
0093   std::cout << "Warning, no recognized jet clustering algorithm provided in FastJetAlgo" << std::endl
0094             << "defaulting to antikt_algorithm" << std::endl;
0095   // return a dummy definition
0096   return fastjet::JetDefinition(fastjet::antikt_algorithm, m_opt.jet_R, fastjet::E_scheme, fastjet::Best);
0097 }
0098 
0099 fastjet::Selector FastJetAlgo::get_selector() const
0100 {
0101   // only selectors available are jet_min_pt and jet_max_eta
0102   if (m_opt.use_jet_max_eta && m_opt.use_jet_min_pt)
0103   {
0104     return fastjet::SelectorAbsRapMax(m_opt.jet_max_eta) && fastjet::SelectorPtMin(m_opt.jet_min_pt);
0105   }
0106   if (m_opt.use_jet_max_eta)
0107   {
0108     return fastjet::SelectorAbsRapMax(m_opt.jet_max_eta);
0109   }
0110 
0111   return fastjet::SelectorPtMin(m_opt.jet_min_pt);
0112 }
0113 
0114 std::vector<fastjet::PseudoJet> FastJetAlgo::cluster_jets(
0115     std::vector<fastjet::PseudoJet>& pseudojets)
0116 {
0117   auto jetdef = get_fastjet_definition();
0118   m_cluseq = new fastjet::ClusterSequence(pseudojets, jetdef);
0119 
0120   if (m_opt.use_jet_selection)
0121   {
0122     auto selector = get_selector();
0123     return fastjet::sorted_by_pt(selector(m_cluseq->inclusive_jets()));
0124   }
0125 
0126   return fastjet::sorted_by_pt(m_cluseq->inclusive_jets());
0127 }
0128 
0129 std::vector<fastjet::PseudoJet> FastJetAlgo::cluster_area_jets(
0130     std::vector<fastjet::PseudoJet>& pseudojets)
0131 {
0132   auto jetdef = get_fastjet_definition();
0133 
0134   fastjet::AreaDefinition area_def(
0135       fastjet::active_area_explicit_ghosts,
0136       fastjet::GhostedAreaSpec(m_opt.ghost_max_rap, 1, m_opt.ghost_area));
0137 
0138   m_cluseqarea = new fastjet::ClusterSequenceArea(pseudojets, jetdef, area_def);
0139 
0140   fastjet::Selector selector = (m_opt.use_jet_selection
0141                                     ? (!fastjet::SelectorIsPureGhost() && get_selector())
0142                                     : !fastjet::SelectorIsPureGhost());
0143 
0144   return fastjet::sorted_by_pt(selector(m_cluseqarea->inclusive_jets()));
0145 }
0146 
0147 float FastJetAlgo::calc_rhomeddens(std::vector<fastjet::PseudoJet>& constituents) const
0148 {
0149   fastjet::AreaDefinition area_def(
0150       fastjet::active_area_explicit_ghosts,
0151       fastjet::GhostedAreaSpec(m_opt.ghost_max_rap, 1, m_opt.ghost_area));
0152 
0153   fastjet::Selector rho_select = (!fastjet::SelectorNHardest(m_opt.nhardestcut_jetmedbkgdens)) * fastjet::SelectorAbsEtaMax(m_opt.etahardestcut_jetmedbkgdens);  // <--
0154 
0155   fastjet::JetDefinition jet_def_bkgd(fastjet::kt_algorithm, m_opt.jet_R);  // <--
0156   fastjet::JetMedianBackgroundEstimator bge{rho_select, jet_def_bkgd, area_def};
0157   bge.set_particles(constituents);
0158   return bge.rho();
0159 }
0160 
0161 std::vector<fastjet::PseudoJet>
0162 FastJetAlgo::jets_to_pseudojets(std::vector<Jet*>& particles) const
0163 {
0164   std::vector<fastjet::PseudoJet> pseudojets;
0165   for (unsigned int ipart = 0; ipart < particles.size(); ++ipart)
0166   {
0167     // fastjet performs strangely with exactly (px,py,pz,E) =
0168     // (0,0,0,0) inputs, such as placeholder towers or those with
0169     // zero'd out energy after CS. this catch also in FastJetAlgoSub
0170 
0171     // Ignore particles with negative/small energies
0172 
0173     if (particles[ipart]->get_e() < m_opt.constituent_min_E)
0174     {
0175       continue;
0176     }
0177     if (!std::isfinite(particles[ipart]->get_px()) ||
0178         !std::isfinite(particles[ipart]->get_py()) ||
0179         !std::isfinite(particles[ipart]->get_pz()) ||
0180         !std::isfinite(particles[ipart]->get_e()))
0181     {
0182       std::cout << PHWHERE << " invalid particle kinematics:"
0183                 << " px: " << particles[ipart]->get_px()
0184                 << " py: " << particles[ipart]->get_py()
0185                 << " pz: " << particles[ipart]->get_pz()
0186                 << " e: " << particles[ipart]->get_e() << std::endl;
0187       gSystem->Exit(1);
0188     }
0189     fastjet::PseudoJet pseudojet(particles[ipart]->get_px(),
0190                                  particles[ipart]->get_py(),
0191                                  particles[ipart]->get_pz(),
0192                                  particles[ipart]->get_e());
0193     if (m_opt.use_constituent_min_pt && pseudojet.perp() < m_opt.constituent_min_pt)
0194     {
0195       continue;
0196     }
0197     pseudojet.set_user_index(ipart);
0198     pseudojets.push_back(pseudojet);
0199   }
0200   return pseudojets;
0201 }
0202 
0203 void FastJetAlgo::first_call_init(JetContainer* jetcont)
0204 {
0205   m_first_cluster_call = false;
0206   m_opt.initialize();
0207 
0208   if (jetcont == nullptr)
0209   {
0210     return;
0211   }
0212 
0213   if (m_opt.doSoftDrop)
0214   {
0215     jetcont->add_property({Jet::PROPERTY::prop_zg, Jet::PROPERTY::prop_Rg, Jet::PROPERTY::prop_mu});
0216     m_zg_index = jetcont->property_index(Jet::PROPERTY::prop_zg);
0217     m_Rg_index = jetcont->property_index(Jet::PROPERTY::prop_Rg);
0218     m_mu_index = jetcont->property_index(Jet::PROPERTY::prop_mu);
0219   }
0220 
0221   if (m_opt.calc_area)
0222   {
0223     jetcont->add_property(Jet::PROPERTY::prop_area);
0224     m_area_index = jetcont->property_index(Jet::PROPERTY::prop_area);
0225   }
0226 
0227   if (m_opt.cs_calc_constsub)
0228   {
0229     cs_bge_rho = new fastjet::GridMedianBackgroundEstimator(m_opt.cs_max_eta, 0.5);
0230     cs_subtractor = new fastjet::contrib::ConstituentSubtractor();
0231     cs_subtractor->set_distance_type(fastjet::contrib::ConstituentSubtractor::deltaR);
0232     cs_subtractor->set_max_distance(m_opt.cs_max_dist);   // free parameter for the maximal allowed distance between particle i and ghost k
0233     cs_subtractor->set_alpha(m_opt.cs_alpha);             // free parameter for the distance measure (the exponent of particle pt). The larger the parameter alpha, the more are favoured the lower pt particles in the subtraction process
0234     cs_subtractor->set_ghost_area(m_opt.cs_ghost_area);   // free parameter for the density of ghosts.
0235     cs_subtractor->set_max_eta(m_opt.cs_max_eta);         // parameter for the maximal eta cut
0236     cs_subtractor->set_background_estimator(cs_bge_rho);  // specify the background estimator to estimate rho.
0237     if (m_opt.cs_max_pt > 0)
0238     {
0239       cs_sel_max_pt = new fastjet::Selector(fastjet::SelectorPtMax(m_opt.cs_max_pt));
0240       cs_subtractor->set_particle_selector(cs_sel_max_pt);
0241     }
0242     cs_subtractor->initialize();
0243     if (m_opt.verbosity > 1)
0244     {
0245       std::cout << cs_subtractor->description() << std::endl;
0246     }
0247   }
0248 
0249   jetcont->set_algo(m_opt.algo);
0250   jetcont->set_jetpar_R(m_opt.jet_R);
0251 }
0252 
0253 void FastJetAlgo::cluster_and_fill(std::vector<Jet*>& particles, JetContainer* jetcont)
0254 {
0255   if (m_first_cluster_call)
0256   {
0257     first_call_init(jetcont);
0258   }
0259   // initalize the properties in JetContainer
0260 
0261   if (m_opt.verbosity > 1)
0262   {
0263     std::cout << "   Verbosity>1 FastJetAlgo::process_event -- entered" << std::endl;
0264   }
0265   if (m_opt.verbosity > 8)
0266   {
0267     std::cout << "   Verbosity>8 #input particles: " << particles.size() << std::endl;
0268   }
0269 
0270   // translate input jets to input fastjets
0271   auto pseudojets = jets_to_pseudojets(particles);
0272 
0273   // if using constituent subtraction, oberve maximum eta and subtract the constituents
0274   if (m_opt.cs_calc_constsub)
0275   {
0276     if (m_opt.verbosity > 100)
0277     {
0278       std::cout << " Before Constituent Subtraction: " << std::endl;
0279       int i = 0;
0280       double sumpt = 0.;
0281       for (const auto& c : pseudojets)
0282       {
0283         sumpt += c.perp();
0284         if (i < 100)
0285         {
0286           std::cout << (boost::format(" jet[%2i] %8.4f  sum %8.4f") % i % c.perp() % sumpt).str() << std::endl;
0287         }
0288         i++;
0289       }
0290       auto _c = pseudojets.back();
0291       std::cout << (boost::format(" jet[%2i] %8.4f  sum %8.4f") % i++ % _c.perp() % sumpt).str() << std::endl
0292                 << std::endl;
0293     }
0294 
0295     pseudojets = fastjet::SelectorAbsEtaMax(m_opt.cs_max_eta)(pseudojets);
0296     cs_bge_rho->set_particles(pseudojets);
0297     auto subtracted_pseudojets = cs_subtractor->subtract_event(pseudojets);
0298     pseudojets = std::move(subtracted_pseudojets);
0299 
0300     if (m_opt.verbosity > 100)
0301     {
0302       std::cout << " After Constituent Subtraction: " << std::endl;
0303       int i = 0;
0304       double sumpt = 0.;
0305       for (const auto& c : pseudojets)
0306       {
0307         sumpt += c.perp();
0308         if (i < 100)
0309         {
0310           std::cout << (boost::format(" jet[%2i] %8.4f  sum %8.4f") % i % c.perp() % sumpt).str() << std::endl;
0311         }
0312         i++;
0313       }
0314       auto _c = pseudojets.back();
0315       std::cout << (boost::format(" jet[%2i] %8.4f  sum %8.4f") % i++ % _c.perp() % sumpt).str() << std::endl
0316                 << std::endl;
0317     }
0318   }
0319 
0320   if (m_opt.calc_jetmedbkgdens)
0321   {
0322     jetcont->set_rho_median(calc_rhomeddens(pseudojets));
0323   }
0324 
0325   auto fastjets = (m_opt.calc_area ? cluster_area_jets(pseudojets) : cluster_jets(pseudojets));
0326 
0327   if (m_opt.verbosity > 8)
0328   {
0329     std::cout << "   Verbosity>8 fastjets: " << fastjets.size() << std::endl;
0330   }
0331   for (unsigned int ijet = 0; ijet < fastjets.size(); ++ijet)
0332   {
0333     auto* jet = jetcont->add_jet();  // put a new Jetv2 into the TClonesArray
0334     jet->set_px(fastjets[ijet].px());
0335     jet->set_py(fastjets[ijet].py());
0336     jet->set_pz(fastjets[ijet].pz());
0337     jet->set_e(fastjets[ijet].e());
0338     jet->set_id(ijet);
0339 
0340     if (m_opt.calc_area)
0341     {
0342       jet->set_property(m_area_index, fastjets[ijet].area());
0343     }
0344 
0345     // if SoftDrop enabled, and jets have > 5 GeV (do not waste time
0346     // on very low-pT jets), run SD and pack output into jet properties
0347     // remove jets that have negative energies
0348     if (m_opt.doSoftDrop && fastjets[ijet].perp() > 5)
0349     {
0350       fastjet::contrib::SoftDrop sd(m_opt.SD_beta, m_opt.SD_zcut);
0351       if (m_opt.verbosity > 5)
0352       {
0353         std::cout << "FastJetAlgo::get_jets : created SoftDrop groomer configuration : "
0354                   << sd.description() << std::endl;
0355       }
0356 
0357       fastjet::PseudoJet sd_jet = sd(fastjets[ijet]);
0358 
0359       if (m_opt.verbosity > 5)
0360       {
0361         std::cout << "original    jet: pt / eta / phi / m = " << fastjets[ijet].perp()
0362                   << " / " << fastjets[ijet].eta() << " / " << fastjets[ijet].phi() << " / "
0363                   << fastjets[ijet].m() << std::endl;
0364         std::cout << "SoftDropped jet: pt / eta / phi / m = " << sd_jet.perp() << " / "
0365                   << sd_jet.eta() << " / " << sd_jet.phi() << " / " << sd_jet.m() << std::endl;
0366 
0367         std::cout << "  delta_R between subjets: " << sd_jet.structure_of<fastjet::contrib::SoftDrop>().delta_R() << std::endl;
0368         std::cout << "  symmetry measure(z):     " << sd_jet.structure_of<fastjet::contrib::SoftDrop>().symmetry() << std::endl;
0369         std::cout << "  mass drop(mu):           " << sd_jet.structure_of<fastjet::contrib::SoftDrop>().mu() << std::endl;
0370       }
0371 
0372       // attach SoftDrop quantities as jet properties
0373       jet->set_property(m_zg_index, sd_jet.structure_of<fastjet::contrib::SoftDrop>().symmetry());
0374       jet->set_property(m_Rg_index, sd_jet.structure_of<fastjet::contrib::SoftDrop>().delta_R());
0375       jet->set_property(m_mu_index, sd_jet.structure_of<fastjet::contrib::SoftDrop>().mu());
0376     }
0377 
0378     // Count clustered components. If desired, put original components into the output jet.
0379     //    int n_clustered = 0;
0380     std::vector<fastjet::PseudoJet> constituents = fastjets[ijet].constituents();
0381     if (m_opt.calc_area)
0382     {
0383       for (auto& comp : constituents)
0384       {
0385         if (comp.is_pure_ghost())
0386         {
0387           continue;
0388         }
0389         //        ++n_clustered;
0390         if (m_opt.save_jet_components)
0391         {
0392           jet->insert_comp(particles[comp.user_index()]->get_comp_vec(), true);
0393         }
0394       }  // end loop over all constituents
0395     }
0396     else
0397     {  // didn't calculate jet area
0398        //      n_clustered += constituents.size();
0399       if (m_opt.save_jet_components)
0400       {
0401         for (auto& comp : constituents)
0402         {
0403           jet->insert_comp(particles[comp.user_index()]->get_comp_vec(), true);
0404         }
0405       }
0406     }
0407     jet->set_comp_sort_flag();  // make surce comp knows it might not be sorted
0408   }
0409   if (m_opt.verbosity > 1)
0410   {
0411     std::cout << "FastJetAlgo::process_event -- exited" << std::endl;
0412   }
0413   delete (m_opt.calc_area ? m_cluseqarea : m_cluseq);  // if (m_cluseq) delete m_cluseq;
0414 }
0415 
0416 std::vector<Jet*> FastJetAlgo::get_jets(std::vector<Jet*> particles)
0417 {
0418   // translate to fastjet
0419   auto pseudojets = jets_to_pseudojets(particles);
0420   auto fastjets = cluster_jets(pseudojets);
0421 
0422   fastjet::contrib::SoftDrop sd(m_opt.SD_beta, m_opt.SD_zcut);
0423   if (m_opt.verbosity > 5)
0424   {
0425     std::cout << "FastJetAlgo::get_jets : created SoftDrop groomer configuration : " << sd.description() << std::endl;
0426   }
0427 
0428   // translate into jet output...
0429   std::vector<Jet*> jets;
0430   for (unsigned int ijet = 0; ijet < fastjets.size(); ++ijet)
0431   {
0432     Jet* jet = new Jetv2();
0433     jet->set_px(fastjets[ijet].px());
0434     jet->set_py(fastjets[ijet].py());
0435     jet->set_pz(fastjets[ijet].pz());
0436     jet->set_e(fastjets[ijet].e());
0437     jet->set_id(ijet);
0438 
0439     // if SoftDrop enabled, and jets have > 5 GeV (do not waste time
0440     // on very low-pT jets), run SD and pack output into jet properties
0441     if (m_opt.doSoftDrop && fastjets[ijet].perp() > 5)
0442     {
0443       fastjet::PseudoJet sd_jet = sd(fastjets[ijet]);
0444 
0445       if (m_opt.verbosity > 5)
0446       {
0447         std::cout << "original    jet: pt / eta / phi / m = " << fastjets[ijet].perp() << " / " << fastjets[ijet].eta() << " / " << fastjets[ijet].phi() << " / " << fastjets[ijet].m() << std::endl;
0448         std::cout << "SoftDropped jet: pt / eta / phi / m = " << sd_jet.perp() << " / " << sd_jet.eta() << " / " << sd_jet.phi() << " / " << sd_jet.m() << std::endl;
0449 
0450         std::cout << "  delta_R between subjets: " << sd_jet.structure_of<fastjet::contrib::SoftDrop>().delta_R() << std::endl;
0451         std::cout << "  symmetry measure(z):     " << sd_jet.structure_of<fastjet::contrib::SoftDrop>().symmetry() << std::endl;
0452         std::cout << "  mass drop(mu):           " << sd_jet.structure_of<fastjet::contrib::SoftDrop>().mu() << std::endl;
0453       }
0454 
0455       // attach SoftDrop quantities as jet properties
0456       jet->set_property(Jet::PROPERTY::prop_zg, sd_jet.structure_of<fastjet::contrib::SoftDrop>().symmetry());
0457       jet->set_property(Jet::PROPERTY::prop_Rg, sd_jet.structure_of<fastjet::contrib::SoftDrop>().delta_R());
0458       jet->set_property(Jet::PROPERTY::prop_mu, sd_jet.structure_of<fastjet::contrib::SoftDrop>().mu());
0459     }
0460 
0461     // Count clustered components. If desired, put original components into the output jet.
0462     //    int n_clustered = 0;
0463     std::vector<fastjet::PseudoJet> constituents = fastjets[ijet].constituents();
0464     if (m_opt.calc_area)
0465     {
0466       for (auto& comp : constituents)
0467       {
0468         if (comp.is_pure_ghost())
0469         {
0470           continue;
0471         }
0472         //        ++n_clustered;
0473         if (m_opt.save_jet_components)
0474         {
0475           jet->insert_comp(particles[comp.user_index()]->get_comp_vec(), true);
0476         }
0477       }  // end loop over all constituents
0478     }
0479     else
0480     {  // didn't save jet area
0481        //      n_clustered += constituents.size();
0482       if (m_opt.save_jet_components)
0483       {
0484         for (auto& comp : constituents)
0485         {
0486           jet->insert_comp(particles[comp.user_index()]->get_comp_vec(), true);
0487         }
0488       }
0489     }
0490     /* jet->set_n_clustered(n_clustered); */
0491     jet->set_comp_sort_flag();  // make surce comp knows it might not be sorted
0492                                 // can alternatively just remove the `true` parameter
0493                                 // from the insert_comp function calls above
0494     jets.push_back(jet);
0495   }
0496 
0497   /* std::sort(jets.begin(), jets.end(), [](Jet*a, Jet*b){ return a->get_et() > b->get_et();} ); */
0498 
0499   if (m_opt.verbosity > 1)
0500   {
0501     std::cout << "FastJetAlgo::process_event -- exited" << std::endl;
0502   }
0503 
0504   delete m_cluseq;
0505 
0506   return jets;
0507 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!