Herwig/HepMCTrigger/HepMCJetTrigger.cc

0001 #include "HepMCJetTrigger.h"
0002
0003 #include <fun4all/SubsysReco.h>
0004 #include <phhepmc/PHHepMCGenEvent.h>
0005 #include <phhepmc/PHHepMCGenEventMap.h>
0006
0007 #include <fun4all/Fun4AllReturnCodes.h>
0008
0009 #include <phool/PHCompositeNode.h>
0010 #include <phool/getClass.h>
0011
0012 #include <fastjet/JetDefinition.hh>
0013
0014 #include <HepMC/GenEvent.h>
0015
0016 #include <fastjet/PseudoJet.hh>
0017 #include <string>
0018 #include <vector>
0019
0020 //____________________________________________________________________________..
0021 // NOLINTNEXTLINE(bugprone-easily-swappable-parameters)
0022 HepMCJetTrigger::HepMCJetTrigger(float trigger_thresh, int n_incom, bool up_lim, const std::string& name)
0023   : SubsysReco(name)
0024   , threshold(trigger_thresh)
0025   , goal_event_number(n_incom)
0026   , set_event_limit(up_lim)
0027 {
0028 }
0029
0030 //____________________________________________________________________________..
0031 int HepMCJetTrigger::process_event(PHCompositeNode* topNode)
0032 {
0033   // std::cout << "HepMCJetTrigger::process_event(PHCompositeNode *topNode) Processing Event" << std::endl;
0034   n_evts++;
0035   if (this->set_event_limit == true)
0036   {  // needed to keep all HepMC output at the same number of events
0037     if (n_good >= this->goal_event_number)
0038     {
0039       return Fun4AllReturnCodes::ABORTEVENT;
0040     }
0041   }
0042   PHHepMCGenEventMap* phg = findNode::getClass<PHHepMCGenEventMap>(topNode, "PHHepMCGenEventMap");
0043   if (!phg)
0044   {
0045     return Fun4AllReturnCodes::ABORTEVENT;
0046   }
0047   for (PHHepMCGenEventMap::ConstIter eventIter = phg->begin(); eventIter != phg->end(); ++eventIter)
0048   {
0049     PHHepMCGenEvent* hepev = eventIter->second;
0050     if (!hepev)
0051     {
0052       return Fun4AllReturnCodes::ABORTEVENT;
0053     }
0054     HepMC::GenEvent* ev = hepev->getEvent();
0055     if (!ev)
0056     {
0057       return Fun4AllReturnCodes::ABORTEVENT;
0058     }
0059     bool const good_event = isGoodEvent(ev);
0060     if (good_event)
0061     {
0062       n_good++;
0063     }
0064     if (!good_event)
0065     {
0066       return Fun4AllReturnCodes::ABORTEVENT;
0067     }
0068   }
0069   return Fun4AllReturnCodes::EVENT_OK;
0070 }
0071
0072 bool HepMCJetTrigger::isGoodEvent(HepMC::GenEvent* e1)
0073 {
0074   // this is really just the call to actually evaluate and return the filter
0075   if (this->threshold == 0)
0076   {
0077     return true;
0078   }
0079   std::vector<fastjet::PseudoJet> const jets = findAllJets(e1);
0080   int const njetsabove = jetsAboveThreshold(jets);
0081   if (njetsabove > 0)
0082   {
0083     return true;
0084   }
0085   return false;
0086 }
0087
0088 std::vector<fastjet::PseudoJet> HepMCJetTrigger::findAllJets(HepMC::GenEvent* e1)
0089 {
0090   // do the fast jet clustering, antikt r=-0.4
0091   fastjet::JetDefinition const jetdef(fastjet::antikt_algorithm, 0.4);
0092   std::vector<fastjet::PseudoJet> input;
0093   std::vector<fastjet::PseudoJet> output;
0094   for (HepMC::GenEvent::particle_const_iterator iter = e1->particles_begin(); iter != e1->particles_end(); ++iter)
0095   {
0096     if (!(*iter)->end_vertex() && (*iter)->status() == 1)
0097     {
0098       auto p = (*iter)->momentum();
0099       fastjet::PseudoJet pj(p.px(), p.py(), p.pz(), p.e());
0100       pj.set_user_index((*iter)->barcode());
0101       input.push_back(pj);
0102     }
0103   }
0104   if (input.empty())
0105   {
0106     return input;
0107   }
0108   fastjet::ClusterSequence const js(input, jetdef);
0109   auto j = js.inclusive_jets();
0110   output.reserve(j.size());
0111   for (const auto& j1 : j)
0112   {
0113     output.push_back(j1);  // just keep in the corect format
0114   }
0115   return output;
0116 }
0117
0118 int HepMCJetTrigger::jetsAboveThreshold(const std::vector<fastjet::PseudoJet>& jets) const
0119 {
0120   // search through for the number of identified jets above the threshold
0121   int n_good_jets = 0;
0122   for (const auto& j : jets)
0123   {
0124     float const pt = j.pt();
0125     if (pt > this->threshold)
0126     {
0127       n_good_jets++;
0128     }
0129   }
0130   return n_good_jets;
0131 }