generators/PHPythia8/PHPy8JetTrigger.cc

0001 #include "PHPy8JetTrigger.h"
0002
0003 #include <Pythia8/Event.h>  // for Event, Particle
0004 #include <Pythia8/Pythia.h>
0005
0006 // fastjet includes
0007 #include <fastjet/ClusterSequence.hh>
0008 #include <fastjet/JetDefinition.hh>
0009 #include <fastjet/PseudoJet.hh>
0010
0011 #include <algorithm>  // for max
0012 #include <cmath>      // for sqrt
0013 #include <cstdlib>    // for abs
0014 #include <iostream>   // for operator<<, endl, basic_ostream
0015 #include <memory>     // for allocator_traits<>::value_type
0016 #include <utility>    // for swap
0017 #include <vector>     // for vector
0018
0019 //__________________________________________________________
0020 PHPy8JetTrigger::PHPy8JetTrigger(const std::string &name)
0021   : PHPy8GenTrigger(name)
0022   , _theEtaHigh(4.0)
0023   , _theEtaLow(1.0)
0024   , _minPt(10.0)
0025   , _minZ(0.0)
0026   , _R(1.0)
0027   , _nconst(0)
0028 {
0029 }
0030
0031 PHPy8JetTrigger::~PHPy8JetTrigger()
0032 {
0033   if (Verbosity() > 0)
0034   {
0035     PrintConfig();
0036   }
0037 }
0038
0039 bool PHPy8JetTrigger::Apply(Pythia8::Pythia *pythia)
0040 {
0041   if (Verbosity() > 2)
0042   {
0043     std::cout << "PHPy8JetTrigger::Apply - pythia event size: "
0044          << pythia->event.size() << std::endl;
0045   }
0046
0047   // Loop over all particles in the event
0048   std::vector<fastjet::PseudoJet> pseudojets;
0049   for (int i = 0; i < pythia->event.size(); ++i)
0050   {
0051     if (pythia->event[i].status() > 0)
0052     {  // only stable particles
0053
0054       // remove some particles (muons, taus, neutrinos)...
0055       // 12 == nu_e
0056       // 13 == muons
0057       // 14 == nu_mu
0058       // 15 == taus
0059       // 16 == nu_tau
0060       if ((abs(pythia->event[i].id()) >= 12) && (abs(pythia->event[i].id()) <= 16))
0061       {
0062         continue;
0063       }
0064
0065       // remove acceptance... _etamin,_etamax
0066       if ((pythia->event[i].px() == 0.0) && (pythia->event[i].py() == 0.0))
0067       {
0068         continue;  // avoid pt=0
0069       }
0070       if ((pythia->event[i].eta() < _theEtaLow ||
0071            pythia->event[i].eta() > _theEtaHigh))
0072       {
0073         continue;
0074       }
0075
0076       fastjet::PseudoJet pseudojet(pythia->event[i].px(),
0077                                    pythia->event[i].py(),
0078                                    pythia->event[i].pz(),
0079                                    pythia->event[i].e());
0080       pseudojet.set_user_index(i);
0081       pseudojets.push_back(pseudojet);
0082     }
0083   }
0084
0085   // Call FastJet
0086
0087   fastjet::JetDefinition *jetdef = new fastjet::JetDefinition(fastjet::antikt_algorithm, _R, fastjet::E_scheme, fastjet::Best);
0088   fastjet::ClusterSequence jetFinder(pseudojets, *jetdef);
0089   std::vector<fastjet::PseudoJet> fastjets = jetFinder.inclusive_jets();
0090   delete jetdef;
0091
0092   bool jetFound = false;
0093   double max_pt = -1;
0094   for (auto &fastjet : fastjets)
0095   {
0096     const double pt = sqrt(fastjet.px() * fastjet.px() + fastjet.py() * fastjet.py());
0097
0098     if (pt > max_pt)
0099     {
0100       max_pt = pt;
0101     }
0102
0103     std::vector<fastjet::PseudoJet> constituents = fastjet.constituents();
0104     int ijet_nconst = constituents.size();
0105
0106     if (pt > _minPt && ijet_nconst >= _nconst)
0107     {
0108       if (_minZ > 0.0)
0109       {
0110         // Loop over constituents, calculate the z of the leading particle
0111
0112         double leading_Z = 0.0;
0113
0114         double jet_ptot = sqrt(fastjet.px() * fastjet.px() +
0115                                fastjet.py() * fastjet.py() +
0116                                fastjet.pz() * fastjet.pz());
0117
0118         for (auto &constituent : constituents)
0119         {
0120           double con_ptot = sqrt(constituent.px() * constituent.px() +
0121                                  constituent.py() * constituent.py() +
0122                                  constituent.pz() * constituent.pz());
0123
0124           double ctheta = (constituent.px() * fastjet.px() +
0125                            constituent.py() * fastjet.py() +
0126                            constituent.pz() * fastjet.pz()) /
0127                           (con_ptot * jet_ptot);
0128
0129           double z_constit = con_ptot * ctheta / jet_ptot;
0130
0131           if (z_constit > leading_Z)
0132           {
0133             leading_Z = z_constit;
0134           }
0135         }
0136
0137         if (leading_Z > _minZ)
0138         {
0139           jetFound = true;
0140           break;
0141         }
0142       }
0143       else
0144       {
0145         jetFound = true;
0146         break;
0147       }
0148     }
0149   }
0150
0151   if (Verbosity() > 2)
0152   {
0153     std::cout << "PHPy8JetTrigger::Apply - max_pt = " << max_pt << ", and jetFound = " << jetFound << std::endl;
0154   }
0155
0156   return jetFound;
0157 }
0158
0159 void PHPy8JetTrigger::SetEtaHighLow(double etaHigh, double etaLow)
0160 {
0161   _theEtaHigh = etaHigh;
0162   _theEtaLow = etaLow;
0163
0164   if (_theEtaHigh < _theEtaLow)
0165   {
0166     std::swap(_theEtaHigh, _theEtaLow);
0167   }
0168 }
0169
0170 void PHPy8JetTrigger::SetMinJetPt(double minPt)
0171 {
0172   _minPt = minPt;
0173 }
0174
0175 void PHPy8JetTrigger::SetMinLeadingZ(double minZ)
0176 {
0177   _minZ = minZ;
0178 }
0179
0180 void PHPy8JetTrigger::SetJetR(double R)
0181 {
0182   _R = R;
0183 }
0184
0185 void PHPy8JetTrigger::SetMinNumConstituents(int nconst)
0186 {
0187   _nconst = nconst;
0188 }
0189
0190 void PHPy8JetTrigger::PrintConfig()
0191 {
0192   std::cout << "---------------- PHPy8JetTrigger::PrintConfig --------------------" << std::endl;
0193
0194   std::cout << "   Particles EtaCut:  " << _theEtaLow << " < eta < " << _theEtaHigh << std::endl;
0195   std::cout << "   Minimum Jet pT: " << _minPt << " GeV/c" << std::endl;
0196   std::cout << "   Anti-kT Radius: " << _R << std::endl;
0197   std::cout << "-----------------------------------------------------------------------" << std::endl;
0198 }