analysis/EICTrigger/PHPy6JetTrigger.C

0001
0002 #include "PHPy6GenTrigger.h"
0003 #include "PHPy6JetTrigger.h"
0004 #include <phool/PHCompositeNode.h>
0005 #include <phool/phool.h>
0006 #include <phool/getClass.h>
0007
0008 #include <phhepmc/PHHepMCGenEvent.h>
0009 #include <HepMC/GenEvent.h>
0010
0011 // fastjet includes
0012 #include <fastjet/JetDefinition.hh>
0013 #include <fastjet/PseudoJet.hh>
0014 #include <fastjet/ClusterSequence.hh>
0015 #include <fastjet/SISConePlugin.hh>
0016 #include <cassert>
0017
0018 using namespace std;
0019
0020 //__________________________________________________________
0021 PHPy6JetTrigger::PHPy6JetTrigger(const std::string &name):
0022   PHPy6GenTrigger(name),
0023   _theEtaHigh(4.0),
0024   _theEtaLow(1.0),
0025   _minPt(10.0),
0026   _R(1.0)
0027  {}
0028
0029 PHPy6JetTrigger::~PHPy6JetTrigger() {
0030   if (_verbosity > 0) PrintConfig();
0031 }
0032
0033 bool PHPy6JetTrigger::Apply(const HepMC::GenEvent* evt) {
0034
0035
0036   // Loop over all particles in the event
0037   int idx = 0;
0038   std::vector<fastjet::PseudoJet> pseudojets;
0039   for ( HepMC::GenEvent::particle_const_iterator p
0040       = evt->particles_begin(); p != evt->particles_end(); ++p ){
0041
0042     idx++;
0043     if ( ((*p)->status()!=1) != 0) continue;
0044
0045     // remove some particles (muons, taus, neutrinos)...
0046     // 12 == nu_e
0047     // 13 == muons
0048     // 14 == nu_mu
0049     // 15 == taus
0050     // 16 == nu_tau
0051     if ((abs(((*p)->pdg_id())) >= 12) && (abs(((*p)->pdg_id())) <= 16)) continue;
0052
0053     // acceptance... _etamin,_etamax
0054     if (((*p)->momentum().px() == 0.0) && ((*p)->momentum().py() == 0.0)) continue; // avoid pt=0
0055     if ( (((*p)->momentum().pseudoRapidity()) < _theEtaLow) ||
0056       (((*p)->momentum().pseudoRapidity()) > _theEtaHigh)) continue;
0057
0058
0059     fastjet::PseudoJet pseudojet ((*p)->momentum().px(),
0060                   (*p)->momentum().py(),
0061                   (*p)->momentum().pz(),
0062                   (*p)->momentum().e());
0063     pseudojet.set_user_index(idx);
0064     pseudojets.push_back(pseudojet);
0065
0066   }
0067
0068   // Call FastJet
0069
0070   fastjet::JetDefinition *jetdef = new fastjet::JetDefinition(fastjet::antikt_algorithm,_R, fastjet::E_scheme,fastjet::Best);
0071   fastjet::ClusterSequence jetFinder(pseudojets,*jetdef);
0072   std::vector<fastjet::PseudoJet> fastjets = jetFinder.inclusive_jets();
0073   delete jetdef;
0074
0075   bool jetFound = false;
0076   double max_pt = -1;
0077   for (unsigned int ijet = 0; ijet < fastjets.size(); ++ijet) {
0078
0079       const double pt =  sqrt(pow(fastjets[ijet].px(),2) + pow(fastjets[ijet].py(),2));
0080
0081       if (pt > max_pt) max_pt = pt;
0082
0083     if(pt > _minPt){
0084       jetFound = true;
0085       break;
0086     }
0087   }
0088
0089   if (_verbosity > 2) {
0090     cout << "PHPy6JetTrigger::Apply - max_pt = "<<max_pt<<", and jetFound = "<<jetFound<<endl;
0091   }
0092
0093   return jetFound;
0094 }
0095
0096 void PHPy6JetTrigger::SetEtaHighLow(double etaHigh, double etaLow) {
0097
0098   _theEtaHigh = etaHigh;
0099   _theEtaLow = etaLow;
0100
0101   if (_theEtaHigh<_theEtaLow)
0102     {
0103       swap(_theEtaHigh, _theEtaLow);
0104     }
0105
0106 }
0107
0108 void PHPy6JetTrigger::SetMinJetPt(double minPt) {
0109   _minPt = minPt;
0110 }
0111
0112 void PHPy6JetTrigger::SetJetR(double R) {
0113   _R = R;
0114 }
0115
0116 void PHPy6JetTrigger::PrintConfig() {
0117   cout << "---------------- PHPy6JetTrigger::PrintConfig --------------------" << endl;
0118
0119   cout << "   Particles EtaCut:  " << _theEtaLow << " < eta < " << _theEtaHigh << endl;
0120   cout << "   Minimum Jet pT: " << _minPt << " GeV/c" << endl;
0121   cout << "   Anti-kT Radius: " << _R << endl;
0122   cout << "-----------------------------------------------------------------------" << endl;
0123 }