3 #include <Pythia8/Event.h>
4 #include <Pythia8/Pythia.h>
7 #include <fastjet/ClusterSequence.hh>
8 #include <fastjet/JetDefinition.hh>
9 #include <fastjet/PseudoJet.hh>
42 cout <<
"PHPy8JetTrigger::Apply - pythia event size: "
43 << pythia->event.size() << endl;
47 std::vector<fastjet::PseudoJet> pseudojets;
48 for (
int i = 0; i < pythia->event.size(); ++i)
50 if (pythia->event[i].status() > 0)
59 if ((
abs(pythia->event[i].id()) >= 12) && (
abs(pythia->event[i].id()) <= 16))
continue;
62 if ((pythia->event[i].px() == 0.0) && (pythia->event[i].py() == 0.0))
continue;
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);
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();
82 bool jetFound =
false;
84 for (
unsigned int ijet = 0; ijet < fastjets.size(); ++ijet)
86 const double pt = sqrt(fastjets[ijet].px() * fastjets[ijet].px() + fastjets[ijet].py() * fastjets[ijet].py());
88 if (pt > max_pt) max_pt =
pt;
90 vector<fastjet::PseudoJet> constituents = fastjets[ijet].constituents();
91 int ijet_nconst = constituents.size();
99 double leading_Z = 0.0;
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());
105 for (
unsigned int j = 0; j < constituents.size(); j++)
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());
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);
116 double z_constit = con_ptot * ctheta / jet_ptot;
118 if (z_constit > leading_Z) leading_Z = z_constit;
121 if (leading_Z >
_minZ)
137 cout <<
"PHPy8JetTrigger::Apply - max_pt = " << max_pt <<
", and jetFound = " << jetFound << endl;
176 cout <<
"---------------- PHPy8JetTrigger::PrintConfig --------------------" << endl;
179 cout <<
" Minimum Jet pT: " <<
_minPt <<
" GeV/c" << endl;
180 cout <<
" Anti-kT Radius: " <<
_R << endl;
181 cout <<
"-----------------------------------------------------------------------" << endl;