Herwig/HepMCTrigger/HepMCParticleTrigger.cc

0001 #include "HepMCParticleTrigger.h"
0002
0003 #include <fun4all/SubsysReco.h>
0004 #include <phhepmc/PHHepMCGenEvent.h>
0005 #include <phhepmc/PHHepMCGenEventMap.h>
0006
0007 #include <fun4all/Fun4AllReturnCodes.h>
0008
0009 #include <phool/PHCompositeNode.h>
0010 #include <phool/getClass.h>
0011
0012 #include <fastjet/JetDefinition.hh>
0013
0014 #include <HepMC/GenEvent.h>
0015
0016 #include <fastjet/PseudoJet.hh>
0017 #include <map>
0018 #include <string>
0019 #include <unordered_set>
0020 #include <unordered_map>
0021 #include <vector>
0022 //____________________________________________________________________________..
0023 // NOLINTNEXTLINE(bugprone-easily-swappable-parameters)
0024 HepMCParticleTrigger::HepMCParticleTrigger(float trigger_thresh, int n_incom, bool up_lim, const std::string& name)
0025   : SubsysReco(name)
0026   , threshold(trigger_thresh)
0027   , goal_event_number(n_incom)
0028   , set_event_limit(up_lim)
0029   , _theEtaHigh(1.1)
0030   , _theEtaLow(-1.1)
0031   , _thePtHigh(999.9)
0032   , _thePtLow(0)
0033   , _thePHigh(999.9)
0034   , _thePLow(-999.9)
0035   , _thePzHigh(999.9)
0036   , _thePzLow(-999.9)
0037   ,
0038
0039     _doEtaHighCut(true)
0040   , _doEtaLowCut(true)
0041   , _doBothEtaCut(true)
0042   ,
0043
0044   _doAbsEtaHighCut(false)
0045   , _doAbsEtaLowCut(false)
0046   , _doBothAbsEtaCut(false)
0047   ,
0048
0049   _doPtHighCut(false)
0050   , _doPtLowCut(false)
0051   , _doBothPtCut(false)
0052   ,
0053     _doPHighCut(false)
0054   , _doPLowCut(false)
0055   , _doBothPCut(false)
0056   ,
0057
0058   _doPzHighCut(false)
0059   , _doPzLowCut(false)
0060   , _doBothPzCut(false)
0061 {
0062   if (threshold != 0)
0063   {
0064     _doPtLowCut = true;
0065     _thePtLow = threshold;
0066   }
0067 }
0068
0069 //____________________________________________________________________________..
0070 int HepMCParticleTrigger::process_event(PHCompositeNode* topNode)
0071 {
0072   // std::cout << "HepMCParticleTrigger::process_event(PHCompositeNode *topNode) Processing Event" << std::endl;
0073   if (this->set_event_limit == true)
0074   {  // needed to keep all HepMC output at the same number of events
0075     if (n_good >= this->goal_event_number)
0076     {
0077       return Fun4AllReturnCodes::ABORTEVENT;
0078     }
0079   }
0080   n_evts++;
0081   bool good_event{false};
0082   PHHepMCGenEventMap* phg = findNode::getClass<PHHepMCGenEventMap>(topNode, "PHHepMCGenEventMap");
0083   if (!phg)
0084   {
0085     return Fun4AllReturnCodes::ABORTEVENT;
0086   }
0087   for (PHHepMCGenEventMap::ConstIter eventIter = phg->begin(); eventIter != phg->end(); ++eventIter)
0088   {
0089     PHHepMCGenEvent* hepev = eventIter->second;
0090     if (!hepev)
0091     {
0092       return Fun4AllReturnCodes::ABORTEVENT;
0093     }
0094     HepMC::GenEvent* ev = hepev->getEvent();
0095     if (!ev)
0096     {
0097       return Fun4AllReturnCodes::ABORTEVENT;
0098     }
0099     good_event = isGoodEvent(ev);
0100     if (!good_event)
0101     {
0102       return Fun4AllReturnCodes::ABORTEVENT;
0103     }
0104   }
0105   if (good_event)
0106   {
0107     n_good++;
0108   }
0109   return Fun4AllReturnCodes::EVENT_OK;
0110 }
0111 void HepMCParticleTrigger::AddParticle(int particlePid)
0112 {
0113   _theParticles.push_back(particlePid);
0114   return;
0115 }
0116 void HepMCParticleTrigger::AddParticles(const std::vector<int>& particles)
0117 {
0118   for (auto p : particles)
0119   {
0120     _theParticles.push_back(p);
0121   }
0122   return;
0123 }
0124
0125 void HepMCParticleTrigger::SetPtHigh(double pt)
0126 {
0127   _thePtHigh = pt;
0128   _doPtHighCut = true;
0129   if (_doPtLowCut)
0130   {
0131     _doBothPtCut = true;
0132   }
0133   return;
0134 }
0135 void HepMCParticleTrigger::SetPtLow(double pt)
0136 {
0137   _thePtLow = pt;
0138   _doPtLowCut = true;
0139   if (_doPtHighCut)
0140   {
0141     _doBothPtCut = true;
0142   }
0143   return;
0144 }
0145 void HepMCParticleTrigger::SetPtHighLow(double ptHigh, double ptLow)
0146 {
0147   _thePtHigh = ptHigh;
0148   _doPtHighCut = true;
0149   _thePtLow = ptLow;
0150   _doPtLowCut = true;
0151   _doBothPtCut = true;
0152   return;
0153 }
0154 void HepMCParticleTrigger::SetPHigh(double pt)
0155 {
0156   _thePHigh = pt;
0157   _doPHighCut = true;
0158   if (_doPLowCut)
0159   {
0160     _doBothPCut = true;
0161   }
0162   return;
0163 }
0164 void HepMCParticleTrigger::SetPLow(double pt)
0165 {
0166   _thePLow = pt;
0167   _doPLowCut = true;
0168   if (_doPHighCut)
0169   {
0170     _doBothPCut = true;
0171   }
0172   return;
0173 }
0174 void HepMCParticleTrigger::SetPHighLow(double ptHigh, double ptLow)
0175 {
0176   _thePHigh = ptHigh;
0177   _doPHighCut = true;
0178   _thePLow = ptLow;
0179   _doPLowCut = true;
0180   _doBothPCut = true;
0181   return;
0182 }
0183 void HepMCParticleTrigger::SetPzHigh(double pt)
0184 {
0185   _thePzHigh = pt;
0186   _doPzHighCut = true;
0187   if (_doPzLowCut)
0188   {
0189     _doBothPzCut = true;
0190   }
0191   return;
0192 }
0193 void HepMCParticleTrigger::SetPzLow(double pt)
0194 {
0195   _thePzLow = pt;
0196   _doPzLowCut = true;
0197   if (_doPzHighCut)
0198   {
0199     _doBothPzCut = true;
0200   }
0201   return;
0202 }
0203 void HepMCParticleTrigger::SetPzHighLow(double ptHigh, double ptLow)
0204 {
0205   _thePzHigh = ptHigh;
0206   _doPzHighCut = true;
0207   _thePzLow = ptLow;
0208   _doPzLowCut = true;
0209   _doBothPzCut = true;
0210   return;
0211 }
0212 void HepMCParticleTrigger::SetEtaHigh(double pt)
0213 {
0214   _theEtaHigh = pt;
0215   _doEtaHighCut = true;
0216   if (_doEtaLowCut)
0217   {
0218     _doBothEtaCut = true;
0219   }
0220   return;
0221 }
0222 void HepMCParticleTrigger::SetEtaLow(double pt)
0223 {
0224   _theEtaLow = pt;
0225   _doEtaLowCut = true;
0226   if (_doEtaHighCut)
0227   {
0228     _doBothEtaCut = true;
0229   }
0230   return;
0231 }
0232 void HepMCParticleTrigger::SetEtaHighLow(double ptHigh, double ptLow)
0233 {
0234   _theEtaHigh = ptHigh;
0235   _doEtaHighCut = true;
0236   _theEtaLow = ptLow;
0237   _doEtaLowCut = true;
0238   _doBothEtaCut = true;
0239   return;
0240 }
0241 void HepMCParticleTrigger::SetAbsEtaHigh(double pt)
0242 {
0243   _theEtaHigh = pt;
0244   _doAbsEtaHighCut = true;
0245   if (_doAbsEtaLowCut)
0246   {
0247     _doBothAbsEtaCut = true;
0248   }
0249   return;
0250 }
0251 void HepMCParticleTrigger::SetAbsEtaLow(double pt)
0252 {
0253   _theEtaLow = pt;
0254   _doAbsEtaLowCut = true;
0255   if (_doAbsEtaHighCut)
0256   {
0257     _doBothAbsEtaCut = true;
0258   }
0259   return;
0260 }
0261 void HepMCParticleTrigger::SetAbsEtaHighLow(double ptHigh, double ptLow)
0262 {
0263   _theEtaHigh = ptHigh;
0264   _doAbsEtaHighCut = true;
0265   _theEtaLow = ptLow;
0266   _doAbsEtaLowCut = true;
0267   _doBothAbsEtaCut = true;
0268   return;
0269 }
0270 bool HepMCParticleTrigger::isGoodEvent(HepMC::GenEvent* e1)
0271 {
0272   std::vector<int> n_trigger_particles = getParticles(e1);
0273   for (auto ntp : n_trigger_particles)
0274   {
0275     if (ntp <= 0)
0276     {
0277       return false;  // make sure all particles have at least 1
0278     }
0279   }
0280   return true;
0281 }
0282
0283 std::vector<int> HepMCParticleTrigger::getParticles(HepMC::GenEvent* e1)
0284 {
0285   std::vector<int> n_trigger{};
0286   std::unordered_set<int> particle_pids;
0287   particle_pids.reserve(_theParticles.size());
0288   for (auto it : _theParticles)
0289   {
0290     particle_pids.insert(std::abs(it));
0291   }
0292   std::unordered_map<int, int> particle_types;
0293   particle_types.reserve(particle_pids.size());
0294
0295   for (HepMC::GenEvent::particle_const_iterator iter = e1->particles_begin(); iter != e1->particles_end(); ++iter)
0296   {
0297     const HepMC::GenParticle *g = *iter;
0298     if (m_doStableParticleOnly && (g->end_vertex() || g->status() != 1))
0299     {
0300       continue;
0301     }
0302
0303     int pid = std::abs(g->pdg_id());
0304     auto ipidx = particle_pids.find(pid);
0305     if(ipidx == particle_pids.end())
0306     {
0307       continue;
0308     }
0309
0310     if (m_rejectFromHadronDecay)
0311     {
0312       if(IsFromHadronDecay(g))
0313       {
0314         continue;
0315       }
0316     }
0317
0318     auto p = g->momentum();
0319     float px = p.px();
0320     float py = p.py();
0321     float pz = p.pz();
0322     float p_M = std::sqrt(std::pow(px, 2) + std::pow(py, 2) + std::pow(pz, 2));
0323     float pt = std::sqrt(std::pow(px, 2) + std::pow(py, 2));
0324     double eta = p.eta();
0325
0326     if ((_doEtaHighCut || _doBothEtaCut) && eta > _theEtaHigh)
0327     {
0328       continue;
0329     }
0330     if ((_doEtaLowCut || _doBothEtaCut) && eta < _theEtaLow)
0331     {
0332       continue;
0333     }
0334     if ((_doAbsEtaHighCut || _doBothAbsEtaCut) && std::abs(eta) > _theEtaHigh)
0335     {
0336       continue;
0337     }
0338     if ((_doAbsEtaLowCut || _doBothAbsEtaCut) && std::abs(eta) < _theEtaLow)
0339     {
0340       continue;
0341     }
0342     if ((_doPtHighCut || _doBothPtCut) && pt > _thePtHigh)
0343     {
0344       continue;
0345     }
0346     if ((_doPtLowCut || _doBothPtCut) && pt < _thePtLow)
0347     {
0348       continue;
0349     }
0350     if ((_doPHighCut || _doBothPCut) && p_M > _thePHigh)
0351     {
0352       continue;
0353     }
0354     if ((_doPLowCut || _doBothPCut) && p_M < _thePLow)
0355     {
0356       continue;
0357     }
0358     if ((_doPzHighCut || _doBothPzCut) && pz > _thePzHigh)
0359     {
0360       continue;
0361     }
0362     if ((_doPzLowCut || _doBothPzCut) && pz < _thePzLow)
0363     {
0364       continue;
0365     }
0366
0367     particle_types[pid]++;
0368
0369     if(particle_types.size() == particle_pids.size())
0370     {
0371       break;
0372     }
0373   }
0374
0375   n_trigger.reserve(_theParticles.size());
0376
0377   for (auto it : _theParticles)
0378   {
0379     auto ptid = particle_types.find(std::abs(it));
0380     n_trigger.push_back((ptid != particle_types.end()) ? ptid->second : 0);
0381   }
0382   return n_trigger;
0383 }
0384
0385 int HepMCParticleTrigger::particleAboveThreshold(const std::map<int, int>& n_particles, int trigger_particle)
0386 {
0387   // search through for the number of identified trigger particles passing cuts
0388   auto it = n_particles.find(std::abs(trigger_particle));
0389   if (it != n_particles.end())
0390   {
0391     return it->second;
0392   }
0393   return 0;
0394 }
0395
0396 bool HepMCParticleTrigger::IsFromHadronDecay(const HepMC::GenParticle* gp)
0397 {
0398   if (!gp)
0399   {
0400     return false;
0401   }
0402
0403   const HepMC::GenVertex* vtx = gp->production_vertex();
0404   if (!vtx)
0405   {
0406     return false;
0407   }
0408
0409   for (auto it = vtx->particles_in_const_begin(); it != vtx->particles_in_const_end(); ++it)
0410   {
0411     const HepMC::GenParticle* mom = *it;
0412     if (!mom)
0413     {
0414       continue;
0415     }
0416
0417     if (IsHadronPDG(mom->pdg_id()))
0418     {
0419       return true;
0420     }
0421   }
0422   return false;
0423 }
0424
0425
0426 bool HepMCParticleTrigger::IsHadronPDG(int _pdg)
0427 {
0428   if(IsIonPDG(_pdg))
0429   {
0430     return false;
0431   }
0432
0433   if(std::abs(_pdg) < 100 )
0434   {
0435     return false;
0436   }
0437
0438   return true;
0439 }