packages/jetbackground/FastJetAlgoSub.cc

0001
0002 #include "FastJetAlgoSub.h"
0003
0004 #include <jetbase/Jet.h>
0005 #include <jetbase/JetContainer.h>
0006
0007 // fastjet includes
0008 #include <fastjet/ClusterSequence.hh>
0009 #include <fastjet/JetDefinition.hh>
0010 #include <fastjet/PseudoJet.hh>
0011
0012 // standard includes
0013 #include <iostream>
0014 #include <sstream>  // for basic_ostringstream
0015 #include <vector>
0016
0017 FastJetAlgoSub::FastJetAlgoSub(const FastJetOptions& options)
0018   : m_opt{options}
0019 {
0020   fastjet::ClusterSequence clusseq;
0021   if (m_opt.verbosity > 0)
0022   {
0023     fastjet::ClusterSequence::print_banner();
0024   }
0025   else
0026   {
0027     std::ostringstream nullstream;
0028     fastjet::ClusterSequence::set_fastjet_banner_stream(&nullstream);
0029     fastjet::ClusterSequence::print_banner();
0030     fastjet::ClusterSequence::set_fastjet_banner_stream(&std::cout);
0031   }
0032 }
0033
0034 void FastJetAlgoSub::identify(std::ostream& os)
0035 {
0036   os << "   FastJetAlgoSub: ";
0037   if (m_opt.algo == Jet::ANTIKT)
0038   {
0039     os << "ANTIKT r=" << m_opt.jet_R;
0040   }
0041   else if (m_opt.algo == Jet::KT)
0042   {
0043     os << "KT r=" << m_opt.jet_R;
0044   }
0045   else if (m_opt.algo == Jet::CAMBRIDGE)
0046   {
0047     os << "CAMBRIDGE r=" << m_opt.jet_R;
0048   }
0049   os << std::endl;
0050 }
0051
0052 /* std::vector<Jet*> FastJetAlgoSub::get_jets(std::vector<Jet*> particles) */
0053 /* { }; //  deprecated by iterating from JetMap; now is JetContainer and most code moved into */
0054 //  into cluster_and_fill
0055
0056 void FastJetAlgoSub::cluster_and_fill(std::vector<Jet*>& particles, JetContainer* jetcont)
0057 {
0058   if (m_opt.verbosity > 1)
0059   {
0060     std::cout << "FastJetAlgoSub::process_event -- entered" << std::endl;
0061   }
0062
0063   // translate to fastjet
0064   std::vector<fastjet::PseudoJet> pseudojets;
0065   for (unsigned int ipart = 0; ipart < particles.size(); ++ipart)
0066   {
0067     float this_e = particles[ipart]->get_e();
0068
0069     if (this_e == 0.)
0070     {
0071       continue;
0072     }
0073
0074     float this_px = particles[ipart]->get_px();
0075     float this_py = particles[ipart]->get_py();
0076     float this_pz = particles[ipart]->get_pz();
0077
0078     if (this_e < 0)
0079     {
0080       // make energy = +1 MeV for purposes of clustering
0081       float e_ratio = 0.001 / this_e;
0082
0083       this_e = this_e * e_ratio;
0084       this_px = this_px * e_ratio;
0085       this_py = this_py * e_ratio;
0086       this_pz = this_pz * e_ratio;
0087
0088       if (m_opt.verbosity > 5)
0089       {
0090         std::cout << " FastJetAlgoSub input particle with negative-E, original kinematics px / py / pz / E = ";
0091         std::cout << particles[ipart]->get_px() << " / " << particles[ipart]->get_py() << " / " << particles[ipart]->get_pz() << " / " << particles[ipart]->get_e() << std::endl;
0092         std::cout << " --> entering with modified kinematics px / py / pz / E = " << this_px << " / " << this_py << " / " << this_pz << " / " << this_e << std::endl;
0093       }
0094     }
0095
0096     fastjet::PseudoJet pseudojet(this_px, this_py, this_pz, this_e);
0097
0098     pseudojet.set_user_index(ipart);
0099     pseudojets.push_back(pseudojet);
0100   }
0101
0102   // run fast jet
0103   fastjet::JetDefinition* jetdef = nullptr;
0104   if (m_opt.algo == Jet::ANTIKT)
0105   {
0106     jetdef = new fastjet::JetDefinition(fastjet::antikt_algorithm, m_opt.jet_R, fastjet::E_scheme, fastjet::Best);
0107   }
0108   else if (m_opt.algo == Jet::KT)
0109   {
0110     jetdef = new fastjet::JetDefinition(fastjet::kt_algorithm, m_opt.jet_R, fastjet::E_scheme, fastjet::Best);
0111   }
0112   else if (m_opt.algo == Jet::CAMBRIDGE)
0113   {
0114     jetdef = new fastjet::JetDefinition(fastjet::cambridge_algorithm, m_opt.jet_R, fastjet::E_scheme, fastjet::Best);
0115   }
0116   else
0117   {
0118     return;
0119   }
0120
0121   fastjet::ClusterSequence jetFinder(pseudojets, *jetdef);
0122   std::vector<fastjet::PseudoJet> fastjets = jetFinder.inclusive_jets();
0123   delete jetdef;
0124
0125   // translate into jet output...
0126   std::vector<Jet*> jets;
0127   for (unsigned int ijet = 0; ijet < fastjets.size(); ++ijet)
0128   {
0129     auto* jet = jetcont->add_jet();
0130
0131     if (m_opt.verbosity > 5 && fastjets[ijet].perp() > 15)
0132     {
0133       std::cout << " FastJetAlgoSub : jet # " << ijet << " comes out of clustering with pt / eta / phi = " << fastjets[ijet].perp() << " / " << fastjets[ijet].eta() << " / " << fastjets[ijet].phi();
0134       std::cout << ", px / py / pz / e = " << fastjets[ijet].px() << " / " << fastjets[ijet].py() << " / " << fastjets[ijet].pz() << " / " << fastjets[ijet].e() << std::endl;
0135     }
0136
0137     float total_px = 0;
0138     float total_py = 0;
0139     float total_pz = 0;
0140     float total_e = 0;
0141
0142     // copy components into output jet
0143     std::vector<fastjet::PseudoJet> comps = fastjets[ijet].constituents();
0144     for (auto& comp : comps)
0145     {
0146       Jet* particle = particles[comp.user_index()];
0147
0148       total_px += particle->get_px();
0149       total_py += particle->get_py();
0150       total_pz += particle->get_pz();
0151       total_e += particle->get_e();
0152       jet->insert_comp(particle->get_comp_vec(), true);
0153     }
0154
0155     jet->set_comp_sort_flag();  // make sure jet know comps might not be sorted
0156                                 // alternatively can just ommit the `true`
0157                                 // in insert_comp call above
0158     jet->set_px(total_px);
0159     jet->set_py(total_py);
0160     jet->set_pz(total_pz);
0161     jet->set_e(total_e);
0162     jet->set_id(ijet);
0163
0164     if (m_opt.verbosity > 5 && fastjets[ijet].perp() > 15)
0165     {
0166       std::cout << " FastJetAlgoSub : jet # " << ijet << " after correcting for proper constituent kinematics, pt / eta / phi = " << jet->get_pt() << " / " << jet->get_eta() << " / " << jet->get_phi();
0167       std::cout << ", px / py / pz / e = " << jet->get_px() << " / " << jet->get_py() << " / " << jet->get_pz() << " / " << jet->get_e() << std::endl;
0168     }
0169   }
0170
0171   if (m_opt.verbosity > 1)
0172   {
0173     std::cout << "FastJetAlgoSub::process_event -- exited" << std::endl;
0174   }
0175   return;
0176 }