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