packages/jetbase/JetReco.cc

0001
0002 #include "JetReco.h"
0003
0004 #include "Jet.h"
0005 #include "JetAlgo.h"
0006 #include "JetContainer.h"
0007 #include "JetContainerv1.h"
0008 #include "JetInput.h"
0009 #include "JetMap.h"
0010 #include "JetMapv1.h"
0011
0012 // PHENIX includes
0013 #include <fun4all/Fun4AllReturnCodes.h>
0014 #include <fun4all/SubsysReco.h>  // for SubsysReco
0015
0016 #include <phool/PHCompositeNode.h>
0017 #include <phool/PHIODataNode.h>
0018 #include <phool/PHNode.h>  // for PHNode
0019 #include <phool/PHNodeIterator.h>
0020 #include <phool/PHObject.h>  // for PHObject
0021 #include <phool/PHTypedNodeIterator.h>
0022 #include <phool/getClass.h>
0023 #include <phool/phool.h>  // for PHWHERE
0024
0025 #include <boost/format.hpp>
0026
0027 // standard includes
0028 #include <cstdlib>  // for exit
0029 #include <fstream>
0030 #include <iostream>
0031 #include <memory>  // for allocator_traits<>::value_type
0032 #include <vector>
0033
0034 JetReco::JetReco(const std::string &name, TRANSITION _which)
0035   : SubsysReco(name)
0036   , which_fill{_which}
0037   , use_jetcon{_which == TRANSITION::JET_CONTAINER || _which == TRANSITION::BOTH || _which == TRANSITION::PRETEND_BOTH}
0038   , use_jetmap{_which == TRANSITION::JET_MAP || _which == TRANSITION::BOTH}
0039 {
0040 }
0041 /* JetReco::JetReco(const std::string &name, TRANSITION _which) */
0042 /*   : SubsysReco(name) */
0043 /*   , which_fill{_which} */
0044 /*   , use_jetcon{false} */
0045 /*   , use_jetmap{true} */
0046 /* { */
0047 /* } */
0048
0049 JetReco::~JetReco()
0050 {
0051   for (auto &_input : _inputs)
0052   {
0053     delete _input;
0054   }
0055   _inputs.clear();
0056   for (auto &_algo : _algos)
0057   {
0058     delete _algo;
0059   }
0060   _algos.clear();
0061   _outputs.clear();
0062 }
0063
0064 int JetReco::InitRun(PHCompositeNode *topNode)
0065 {
0066   if (Verbosity() > 0)
0067   {
0068     std::cout << "========================== JetReco::InitRun() =============================" << std::endl;
0069     std::cout << " Input Selections:" << std::endl;
0070     for (auto &_input : _inputs)
0071     {
0072       _input->identify();
0073     }
0074     std::cout << " Algorithms:" << std::endl;
0075     for (auto &_algo : _algos)
0076     {
0077       _algo->identify();
0078     }
0079     std::cout << "===========================================================================" << std::endl;
0080   }
0081
0082   return CreateNodes(topNode);
0083 }
0084
0085 int JetReco::process_event(PHCompositeNode *topNode)
0086 {
0087   if (Verbosity() > 1)
0088   {
0089     std::cout << "JetReco::process_event -- entered" << std::endl;
0090   }
0091
0092   //------------------------------------------------------------------
0093   // This will also need to go into TClonesArrays in a future revision
0094   // Get Objects off of the Node Tree
0095   //------------------------------------------------------------------
0096
0097   std::vector<Jet *> inputs;  // owns memory
0098   for (auto &_input : _inputs)
0099   {
0100     std::vector<Jet *> parts = _input->get_input(topNode);
0101     for (auto &part : parts)
0102     {
0103       inputs.push_back(part);
0104       inputs.back()->set_id(inputs.size() - 1);  // unique ids ensured
0105     }
0106   }
0107
0108   //---------------------------
0109   // Run the jet reconstruction
0110   //---------------------------
0111   for (unsigned int ialgo = 0; ialgo < _algos.size(); ++ialgo)
0112   {
0113     // send the output somewhere on the DST
0114     /* if (_fill_JetContainer) { */
0115     if (use_jetcon)
0116     {
0117       if (Verbosity() > 5)
0118       {
0119         std::cout << " Verbosity>5:: filling JetContainter for " << JC_name(_outputs[ialgo]) << std::endl;
0120       }
0121       FillJetContainer(topNode, ialgo, inputs);
0122     }
0123     if (use_jetmap)
0124     {
0125       if (Verbosity() > 5)
0126       {
0127         std::cout << " Verbosity>5:: filling jetnode for " << _outputs[ialgo] << std::endl;
0128       }
0129       std::vector<Jet *> jets = _algos[ialgo]->get_jets(inputs);  // owns memory
0130       FillJetNode(topNode, ialgo, jets);
0131     }
0132
0133     if (false)
0134     {  // These printouts were used in comparing the two map methods
0135        // It can be removed later
0136       if (use_jetmap)
0137       {
0138         std::fstream fout;
0139         fout.open("fixme_jetmap", std::fstream::app);
0140         fout << " Printing jet results " << std::endl;
0141         JetMap *jetmap = findNode::getClass<JetMap>(topNode, _outputs[ialgo]);
0142         int ifixme = 0;
0143         for (auto &_jet : *jetmap)
0144         {
0145           auto *jet = _jet.second;
0146           fout << (boost::format(" jet[%2i] ncon:pt:eta:phi [%6i,%6.3f,%6.3f,%6.3f]") % (++ifixme) % ((int) jet->size_comp()) % jet->get_pt() % jet->get_eta() % jet->get_phi()).str() << std::endl;
0147           std::vector<std::pair<int, int>> vconst;
0148           /*legacy*/ for (auto _comp = jet->begin_comp(); _comp != jet->end_comp(); ++_comp)
0149           {
0150             vconst.emplace_back((int) _comp->first, (int) _comp->second);
0151           }
0152           std::sort(vconst.begin(), vconst.end(), [](std::pair<int, int> a, std::pair<int, int> b)
0153                     {
0154               if      (a.first == b.first) { return a.second < b.second;
0155               } return a.first < b.first; });
0156           fout << " constituents: ";
0157           int iconst = 0;
0158           for (const auto &p : vconst)
0159           {
0160             fout << (boost::format("  c(%3i) %3i->%3i") % (iconst++) % p.first % p.second).str() << std::endl;
0161           }
0162           fout << std::endl;
0163         }
0164       }
0165       if (use_jetcon)
0166       {
0167         std::fstream fout;
0168         fout.open("fixme_jetcont", std::fstream::app);
0169         fout << " Printing jet results" << std::endl;
0170         JetContainer *jet_cont = findNode::getClass<JetContainer>(topNode, JC_name(_outputs[ialgo]));
0171         int ifixme = 0;
0172
0173         /* for (auto jet = jet_cont->begin(); jet != jet_cont->end(); ++jet) { */
0174         for (auto *jet : *jet_cont)
0175         {
0176           fout << (boost::format(" jet[%2i] ncon:pt:eta:phi [%6i,%6.3f,%6.3f,%6.3f]") % (++ifixme) % ((int) jet->size_comp()) % jet->get_pt() % jet->get_eta() % jet->get_phi()).str() << std::endl;
0177
0178           fout << " constituents: ";
0179           int iconst = 0;
0180           for (const auto &p : jet->get_comp_vec())
0181           {
0182             fout << (boost::format("  c(%3i) %3i->%3i") % (iconst++) % p.first % p.second).str() << std::endl;
0183           }
0184
0185           fout << std::endl;
0186         }
0187       }
0188     }
0189   }
0190
0191   // clean up input vector
0192   // <- another place where TClonesArray's would make this more efficient
0193   for (auto &input : inputs)
0194   {
0195     delete input;
0196   }
0197   inputs.clear();
0198
0199   if (Verbosity() > 1)
0200   {
0201     std::cout << "JetReco::process_event -- exited" << std::endl;
0202   }
0203
0204   return Fun4AllReturnCodes::EVENT_OK;
0205 }
0206
0207 int JetReco::CreateNodes(PHCompositeNode *topNode)
0208 {
0209   PHNodeIterator iter(topNode);
0210
0211   // Looking for the DST node
0212   PHCompositeNode *dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
0213   if (!dstNode)
0214   {
0215     std::cout << PHWHERE << "DST Node missing, doing nothing." << std::endl;
0216     return Fun4AllReturnCodes::ABORTRUN;
0217   }
0218
0219   // Create the AntiKt node if required
0220   PHCompositeNode *AlgoNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", _algonode));
0221   if (!AlgoNode)
0222   {
0223     AlgoNode = new PHCompositeNode(_algonode);
0224     dstNode->addNode(AlgoNode);
0225   }
0226
0227   // Create the Input node if required
0228   PHCompositeNode *InputNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", _inputnode));
0229   if (!InputNode)
0230   {
0231     InputNode = new PHCompositeNode(_inputnode);
0232     AlgoNode->addNode(InputNode);
0233   }
0234
0235   for (auto &_output : _outputs)
0236   {
0237     if (use_jetcon)
0238     {
0239       JetContainer *jetconn = findNode::getClass<JetContainer>(topNode, JC_name(_output));
0240       if (!jetconn)
0241       {
0242         jetconn = new JetContainerv1();
0243         PHIODataNode<PHObject> *JetContainerNode = new PHIODataNode<PHObject>(jetconn, JC_name(_output), "PHObject");
0244         InputNode->addNode(JetContainerNode);
0245       }
0246     }
0247     if (use_jetmap)
0248     {
0249       JetMap *jets = findNode::getClass<JetMap>(topNode, _output);
0250       if (!jets)
0251       {
0252         jets = new JetMapv1();
0253         PHIODataNode<PHObject> *JetMapNode = new PHIODataNode<PHObject>(jets, _output, "PHObject");
0254         InputNode->addNode(JetMapNode);
0255       }
0256     }
0257   }
0258
0259   return Fun4AllReturnCodes::EVENT_OK;
0260 }
0261
0262 void JetReco::FillJetNode(PHCompositeNode *topNode, int ipos, const std::vector<Jet *> &jets)
0263 {
0264   JetMap *jetmap = findNode::getClass<JetMap>(topNode, _outputs[ipos]);
0265   if (!jetmap)
0266   {
0267     std::cout << PHWHERE << " ERROR: Can't find JetMap: " << _outputs[ipos] << std::endl;
0268     exit(-1);
0269   }
0270   jetmap->Reset();
0271   jetmap->set_algo(_algos[ipos]->get_algo());
0272   jetmap->set_par(_algos[ipos]->get_par());
0273   for (auto &_input : _inputs)
0274   {
0275     jetmap->insert_src(_input->get_src());
0276   }
0277
0278   for (const auto &jet : jets)
0279   {
0280     jetmap->insert(jet);  // map takes ownership, sets unique id
0281   }
0282
0283   return;
0284 }
0285
0286 void JetReco::FillJetContainer(PHCompositeNode *topNode, int ipos, std::vector<Jet *> &inputs)
0287 {
0288   JetContainer *jetconn = findNode::getClass<JetContainer>(topNode, JC_name(_outputs[ipos]));
0289   if (!jetconn)
0290   {
0291     std::cout << PHWHERE << " ERROR: Can't find JetContainer: " << _outputs[ipos] << std::endl;
0292     exit(-1);
0293   }
0294   jetconn->Reset();
0295   _algos[ipos]->cluster_and_fill(inputs, jetconn);  // fills the jet container with clustered jets
0296   for (auto &_input : _inputs)
0297   {
0298     jetconn->insert_src(_input->get_src());
0299   }
0300
0301   if (Verbosity() > 7)
0302   {
0303     std::cout << " Verbosity()>7:: jets in container " << _outputs[ipos] << std::endl;
0304     jetconn->print_jets();
0305   }
0306
0307   return;
0308 }
0309
0310 JetAlgo *JetReco::get_algo(unsigned int which_algo)
0311 {
0312   if (_algos.empty())
0313   {
0314     std::cout << PHWHERE << std::endl
0315               << " JetReco has only " << _algos.size() << " JetAlgos; cannot get the one indexed " << which_algo << std::endl;
0316     assert(which_algo < _algos.size());
0317     return nullptr;
0318   }
0319   return _algos[which_algo];
0320 }