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     float w_t_sum = 0;
0142     float w_e_sum = 0;
0143     // copy components into output jet
0144     std::vector<fastjet::PseudoJet> comps = fastjets[ijet].constituents();
0145     for (auto& comp : comps)
0146     {
0147       Jet* particle = particles[comp.user_index()];
0148
0149       total_px += particle->get_px();
0150       total_py += particle->get_py();
0151       total_pz += particle->get_pz();
0152       total_e += particle->get_e();
0153       if(particle->size_properties() > Jet::PROPERTY::prop_t)
0154     {
0155       if(!std::isnan(particle->get_property(Jet::PROPERTY::prop_t)))
0156         {
0157           w_t_sum += particle->get_property(Jet::PROPERTY::prop_t) * particle->get_e();
0158           w_e_sum += particle->get_e();
0159         }
0160     }
0161     jet->insert_comp(particle->get_comp_vec(), true);
0162     }
0163     if(jet->size_properties() < Jet::PROPERTY::prop_t+1)
0164       {
0165     jet->resize_properties(Jet::PROPERTY::prop_t + 1);
0166       }
0167     jet->set_property(Jet::PROPERTY::prop_t,w_t_sum/w_e_sum); //This intentionally becomes nan (0/0) if the
0168                                                               //value has not been filled at all so that
0169                                                               //the jet auto-fails timing cuts if necessary
0170
0171     jet->set_comp_sort_flag();  // make sure jet know comps might not be sorted
0172                                 // alternatively can just ommit the `true`
0173                                 // in insert_comp call above
0174     jet->set_px(total_px);
0175     jet->set_py(total_py);
0176     jet->set_pz(total_pz);
0177     jet->set_e(total_e);
0178     jet->set_id(ijet);
0179
0180     if (m_opt.verbosity > 5 && fastjets[ijet].perp() > 15)
0181     {
0182       std::cout << " FastJetAlgoSub : jet # " << ijet << " after correcting for proper constituent kinematics, pt / eta / phi = " << jet->get_pt() << " / " << jet->get_eta() << " / " << jet->get_phi();
0183       std::cout << ", px / py / pz / e = " << jet->get_px() << " / " << jet->get_py() << " / " << jet->get_pz() << " / " << jet->get_e() << std::endl;
0184     }
0185   }
0186
0187   if (m_opt.verbosity > 1)
0188   {
0189     std::cout << "FastJetAlgoSub::process_event -- exited" << std::endl;
0190   }
0191   return;
0192 }