ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHPy8JetTrigger.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHPy8JetTrigger.cc
1 #include "PHPy8JetTrigger.h"
2 
3 #include <Pythia8/Event.h> // for Event, Particle
4 #include <Pythia8/Pythia.h>
5 
6 // fastjet includes
7 #include <fastjet/ClusterSequence.hh>
8 #include <fastjet/JetDefinition.hh>
9 #include <fastjet/PseudoJet.hh>
10 
11 #include <algorithm> // for max
12 #include <cmath> // for sqrt
13 #include <cstdlib> // for abs
14 #include <iostream> // for operator<<, endl, basic_ostream
15 #include <memory> // for allocator_traits<>::value_type
16 #include <utility> // for swap
17 #include <vector> // for vector
18 
19 using namespace std;
20 
21 //__________________________________________________________
23  : PHPy8GenTrigger(name)
24  , _theEtaHigh(4.0)
25  , _theEtaLow(1.0)
26  , _minPt(10.0)
27  , _minZ(0.0)
28  , _R(1.0)
29  , _nconst(0)
30 {
31 }
32 
34 {
35  if (Verbosity() > 0) PrintConfig();
36 }
37 
38 bool PHPy8JetTrigger::Apply(Pythia8::Pythia *pythia)
39 {
40  if (Verbosity() > 2)
41  {
42  cout << "PHPy8JetTrigger::Apply - pythia event size: "
43  << pythia->event.size() << endl;
44  }
45 
46  // Loop over all particles in the event
47  std::vector<fastjet::PseudoJet> pseudojets;
48  for (int i = 0; i < pythia->event.size(); ++i)
49  {
50  if (pythia->event[i].status() > 0)
51  { //only stable particles
52 
53  // remove some particles (muons, taus, neutrinos)...
54  // 12 == nu_e
55  // 13 == muons
56  // 14 == nu_mu
57  // 15 == taus
58  // 16 == nu_tau
59  if ((abs(pythia->event[i].id()) >= 12) && (abs(pythia->event[i].id()) <= 16)) continue;
60 
61  // remove acceptance... _etamin,_etamax
62  if ((pythia->event[i].px() == 0.0) && (pythia->event[i].py() == 0.0)) continue; // avoid pt=0
63  if ((pythia->event[i].eta() < _theEtaLow ||
64  pythia->event[i].eta() > _theEtaHigh)) continue;
65 
66  fastjet::PseudoJet pseudojet(pythia->event[i].px(),
67  pythia->event[i].py(),
68  pythia->event[i].pz(),
69  pythia->event[i].e());
70  pseudojet.set_user_index(i);
71  pseudojets.push_back(pseudojet);
72  }
73  }
74 
75  // Call FastJet
76 
77  fastjet::JetDefinition *jetdef = new fastjet::JetDefinition(fastjet::antikt_algorithm, _R, fastjet::E_scheme, fastjet::Best);
78  fastjet::ClusterSequence jetFinder(pseudojets, *jetdef);
79  std::vector<fastjet::PseudoJet> fastjets = jetFinder.inclusive_jets();
80  delete jetdef;
81 
82  bool jetFound = false;
83  double max_pt = -1;
84  for (unsigned int ijet = 0; ijet < fastjets.size(); ++ijet)
85  {
86  const double pt = sqrt(fastjets[ijet].px() * fastjets[ijet].px() + fastjets[ijet].py() * fastjets[ijet].py());
87 
88  if (pt > max_pt) max_pt = pt;
89 
90  vector<fastjet::PseudoJet> constituents = fastjets[ijet].constituents();
91  int ijet_nconst = constituents.size();
92 
93  if (pt > _minPt && ijet_nconst >= _nconst)
94  {
95  if (_minZ > 0.0)
96  {
97  // Loop over constituents, calculate the z of the leading particle
98 
99  double leading_Z = 0.0;
100 
101  double jet_ptot = sqrt(fastjets[ijet].px() * fastjets[ijet].px() +
102  fastjets[ijet].py() * fastjets[ijet].py() +
103  fastjets[ijet].pz() * fastjets[ijet].pz());
104 
105  for (unsigned int j = 0; j < constituents.size(); j++)
106  {
107  double con_ptot = sqrt(constituents[j].px() * constituents[j].px() +
108  constituents[j].py() * constituents[j].py() +
109  constituents[j].pz() * constituents[j].pz());
110 
111  double ctheta = (constituents[j].px() * fastjets[ijet].px() +
112  constituents[j].py() * fastjets[ijet].py() +
113  constituents[j].pz() * fastjets[ijet].pz()) /
114  (con_ptot * jet_ptot);
115 
116  double z_constit = con_ptot * ctheta / jet_ptot;
117 
118  if (z_constit > leading_Z) leading_Z = z_constit;
119  }
120 
121  if (leading_Z > _minZ)
122  {
123  jetFound = true;
124  break;
125  }
126  }
127  else
128  {
129  jetFound = true;
130  break;
131  }
132  }
133  }
134 
135  if (Verbosity() > 2)
136  {
137  cout << "PHPy8JetTrigger::Apply - max_pt = " << max_pt << ", and jetFound = " << jetFound << endl;
138  }
139 
140  return jetFound;
141 }
142 
143 void PHPy8JetTrigger::SetEtaHighLow(double etaHigh, double etaLow)
144 {
145  _theEtaHigh = etaHigh;
146  _theEtaLow = etaLow;
147 
148  if (_theEtaHigh < _theEtaLow)
149  {
151  }
152 }
153 
155 {
156  _minPt = minPt;
157 }
158 
160 {
161  _minZ = minZ;
162 }
163 
165 {
166  _R = R;
167 }
168 
170 {
171  _nconst = nconst;
172 }
173 
175 {
176  cout << "---------------- PHPy8JetTrigger::PrintConfig --------------------" << endl;
177 
178  cout << " Particles EtaCut: " << _theEtaLow << " < eta < " << _theEtaHigh << endl;
179  cout << " Minimum Jet pT: " << _minPt << " GeV/c" << endl;
180  cout << " Anti-kT Radius: " << _R << endl;
181  cout << "-----------------------------------------------------------------------" << endl;
182 }