ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHPy6JetTrigger.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHPy6JetTrigger.cc
1 #include "PHPy6JetTrigger.h"
2 #include "PHPy6GenTrigger.h"
3 
4 #include <HepMC/GenEvent.h>
5 #include <HepMC/GenParticle.h> // for GenParticle
6 #include <HepMC/SimpleVector.h> // for FourVector
7 
8 // fastjet includes
9 #include <fastjet/ClusterSequence.hh>
10 #include <fastjet/JetDefinition.hh>
11 #include <fastjet/PseudoJet.hh>
12 
13 #include <cmath> // for sqrt
14 #include <cstdlib> // for abs
15 #include <iostream> // for operator<<, endl, basic_ostream
16 #include <memory> // for allocator_traits<>::value_type
17 #include <utility> // for swap
18 #include <vector> // for vector
19 
20 using namespace std;
21 
22 //__________________________________________________________
24  : PHPy6GenTrigger(name)
25 {
26 }
27 
29 {
30  if (Verbosity() > 0) PrintConfig();
31 }
32 
33 bool PHPy6JetTrigger::Apply(const HepMC::GenEvent *evt)
34 {
35  // Loop over all particles in the event
36  int idx = 0;
37  std::vector<fastjet::PseudoJet> pseudojets;
38  for (HepMC::GenEvent::particle_const_iterator p = evt->particles_begin(); p != evt->particles_end(); ++p)
39  {
40  idx++;
41  if (((*p)->status() != 1) != 0) continue;
42 
43  // remove some particles (muons, taus, neutrinos)...
44  // 12 == nu_e
45  // 13 == muons
46  // 14 == nu_mu
47  // 15 == taus
48  // 16 == nu_tau
49  if ((abs(((*p)->pdg_id())) >= 12) && (abs(((*p)->pdg_id())) <= 16)) continue;
50 
51  // acceptance... _etamin,_etamax
52  if (((*p)->momentum().px() == 0.0) && ((*p)->momentum().py() == 0.0)) continue; // avoid pt=0
53  if ((((*p)->momentum().pseudoRapidity()) < m_theEtaLow) ||
54  (((*p)->momentum().pseudoRapidity()) > m_theEtaHigh)) continue;
55 
56  fastjet::PseudoJet pseudojet((*p)->momentum().px(),
57  (*p)->momentum().py(),
58  (*p)->momentum().pz(),
59  (*p)->momentum().e());
60  pseudojet.set_user_index(idx);
61  pseudojets.push_back(pseudojet);
62  }
63 
64  // Call FastJet
65 
66  fastjet::JetDefinition *jetdef = new fastjet::JetDefinition(fastjet::antikt_algorithm, m_R, fastjet::E_scheme, fastjet::Best);
67  fastjet::ClusterSequence jetFinder(pseudojets, *jetdef);
68  std::vector<fastjet::PseudoJet> fastjets = jetFinder.inclusive_jets();
69  delete jetdef;
70 
71  bool jetFound = false;
72  double max_pt = -1;
73  for (unsigned int ijet = 0; ijet < fastjets.size(); ++ijet)
74  {
75  const double pt = sqrt(pow(fastjets[ijet].px(), 2) + pow(fastjets[ijet].py(), 2));
76 
77  if (pt > max_pt) max_pt = pt;
78 
79  vector<fastjet::PseudoJet> constituents = fastjets[ijet].constituents();
80  const int nconst = constituents.size();
81 
82  if (pt > m_minPt && nconst >= m_nconst)
83  {
84  jetFound = true;
85  break;
86  }
87  }
88 
89  if (Verbosity() > 2)
90  {
91  cout << "PHPy6JetTrigger::Apply - max_pt = " << max_pt << ", and jetFound = " << jetFound << endl;
92  }
93 
94  return jetFound;
95 }
96 
97 void PHPy6JetTrigger::SetEtaHighLow(double etaHigh, double etaLow)
98 {
99  m_theEtaHigh = etaHigh;
100  m_theEtaLow = etaLow;
101 
103  {
105  }
106 }
107 
109 {
110  cout << "---------------- PHPy6JetTrigger::PrintConfig --------------------" << endl;
111 
112  cout << " Particles EtaCut: " << m_theEtaLow << " < eta < " << m_theEtaHigh << endl;
113  cout << " Minimum Jet pT: " << m_minPt << " GeV/c" << endl;
114  cout << " Anti-kT Radius: " << m_R << endl;
115  cout << "-----------------------------------------------------------------------" << endl;
116 }