generators/PHHerwig7/PHHerwig7JetTrigger.cc

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