packages/jetbase/JetProbeMaker.cc

0001 #include "JetProbeMaker.h"
0002
0003 #include "JetContainer.h"
0004 #include "JetContainerv1.h"
0005 #include "Jetv2.h"
0006
0007 #include <fun4all/Fun4AllReturnCodes.h>
0008 #include <phool/PHCompositeNode.h>
0009 #include <phool/PHIODataNode.h>
0010 #include <phool/PHNode.h>  // for PHNode
0011 #include <phool/PHNodeIterator.h>
0012 #include <phool/PHObject.h>      // for PHObject
0013 #include <phool/PHRandomSeed.h>  // for PHRandomSeed
0014 #include <phool/getClass.h>
0015 #include <phool/phool.h>  // for PHWHERE
0016
0017 #include <TRandom3.h>
0018 #include <fastjet/PseudoJet.hh>
0019
0020 int JetProbeMaker::process_event(PHCompositeNode * /*topNode*/)
0021 {
0022   // update the jet probe
0023   float phi = M_PI * (1. - 2. * gsl_rng_uniform(m_rng.get()));
0024   float eta = _eta_min + (gsl_rng_uniform(m_rng.get()) * _eta_range);
0025   float pt = (_pt_range == 0.
0026                   ? _pt_min
0027                   : _pt_min + (gsl_rng_uniform(m_rng.get()) * _pt_range));
0028
0029   fastjet::PseudoJet fjet{};
0030   fjet.reset_PtYPhiM(pt, eta, phi);
0031
0032   Jetv2 *jet = (Jetv2 *) _jets->add_jet();
0033   jet->set_px(fjet.px());
0034   jet->set_py(fjet.py());
0035   jet->set_pz(fjet.pz());
0036   jet->set_e(fjet.e());
0037   jet->insert_comp(Jet::SRC::JET_PROBE, 0, false);
0038
0039   return Fun4AllReturnCodes::EVENT_OK;
0040 }
0041
0042 JetProbeMaker::JetProbeMaker(const std::string &name)
0043   : SubsysReco(name)
0044 {
0045   // initialize rng
0046   const uint seed = PHRandomSeed();
0047   m_rng.reset(gsl_rng_alloc(gsl_rng_mt19937));
0048   gsl_rng_set(m_rng.get(), seed);
0049 };
0050
0051 int JetProbeMaker::InitRun(PHCompositeNode *topNode)
0052 {
0053   // Create the Input node if required
0054   PHNodeIterator iter(topNode);
0055
0056   // Looking for the DST node
0057   PHCompositeNode *dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
0058   if (!dstNode)
0059   {
0060     std::cout << PHWHERE << "DST Node missing, doing nothing." << std::endl;
0061     return Fun4AllReturnCodes::ABORTRUN;
0062   }
0063
0064   // Create the Input node if required
0065   PHCompositeNode *probeNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "JetProbeNode"));
0066   if (!probeNode)
0067   {
0068     probeNode = new PHCompositeNode("JetProbeNode");
0069     dstNode->addNode(probeNode);
0070   }
0071
0072   _jets = findNode::getClass<JetContainer>(topNode, "JetProbeContainer");
0073   if (!_jets)
0074   {
0075     _jets = new JetContainerv1();
0076     PHIODataNode<PHObject> *JetContainerNode = new PHIODataNode<PHObject>(_jets, "JetProbeContainer", "PHObject");
0077     dstNode->addNode(JetContainerNode);
0078   }
0079
0080   return Fun4AllReturnCodes::EVENT_OK;
0081 }