g4simulation/g4jets/TruthJetInput.cc

0001
0002 #include "TruthJetInput.h"
0003
0004 #include <g4main/PHG4Particle.h>
0005 #include <g4main/PHG4TruthInfoContainer.h>
0006
0007 #include <jetbase/Jet.h>
0008 #include <jetbase/Jetv2.h>
0009
0010 #include <phool/getClass.h>
0011 #include <phool/phool.h>  // for PHWHERE
0012
0013 // standard includes
0014 #include <algorithm>  // std::find
0015 #include <cmath>      // for asinh, sqrt
0016 #include <cstdlib>
0017 #include <iostream>
0018 #include <map>      // for _Rb_tree_const_iterator
0019 #include <utility>  // for pair
0020 #include <vector>
0021
0022 TruthJetInput::TruthJetInput(Jet::SRC input)
0023   : m_Input(input)
0024 {
0025   return;
0026 }
0027
0028 void TruthJetInput::identify(std::ostream &os)
0029 {
0030   os << "   TruthJetInput: G4TruthInfo to Jet::PARTICLE";
0031   if (use_embed_stream())
0032   {
0033     os << ". Processing embedded streams: ";
0034     for (int it : m_EmbedID)
0035     {
0036       os << it << ", ";
0037     }
0038   }
0039   os << std::endl;
0040 }
0041
0042 std::vector<Jet *> TruthJetInput::get_input(PHCompositeNode *topNode)
0043 {
0044   if (Verbosity() > 0)
0045   {
0046     std::cout << "TruthJetInput::process_event -- entered" << std::endl;
0047   }
0048
0049   // Pull the reconstructed track information off the node tree...
0050   PHG4TruthInfoContainer *truthinfo = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
0051   if (!truthinfo)
0052   {
0053     std::cout << PHWHERE << " ERROR: Can't find G4TruthInfo" << std::endl;
0054     return std::vector<Jet *>();
0055   }
0056
0057   // FIXME
0058   //   - This requires care: selecting only charged primary particles will
0059   //     miss the decay products of resonances like D0, B0, etc.
0060   //   - For charged particle input, it might be better to
0061   //     (1) Loop on GetPrimaryParticleRange
0062   //     (2) Or add a loop on GetSecondaryParticleRange
0063   std::vector<Jet *> pseudojets;
0064   PHG4TruthInfoContainer::ConstRange range = truthinfo->GetSPHENIXPrimaryParticleRange();
0065   for (PHG4TruthInfoContainer::ConstIterator iter = range.first;
0066        iter != range.second;
0067        ++iter)
0068   {
0069     PHG4Particle *part = iter->second;
0070
0071     if (use_embed_stream())
0072     {
0073       const int this_embed_id = truthinfo->isEmbeded(part->get_track_id());
0074
0075       if (std::find(m_EmbedID.begin(), m_EmbedID.end(), this_embed_id) == m_EmbedID.end())
0076       {
0077         continue;  // reject particle as it is not in the interested embedding stream.
0078       }
0079     }
0080
0081     // remove some particles (muons, taus, neutrinos)...
0082     // 12 == nu_e
0083     // 13 == muons
0084     // 14 == nu_mu
0085     // 15 == taus
0086     // 16 == nu_tau
0087     if ((abs(part->get_pid()) >= 12) && (abs(part->get_pid()) <= 16))
0088     {
0089       continue;
0090     }
0091
0092     // if looking at only charged particles, remove neutrals
0093     if (m_Input == Jet::CHARGED_PARTICLE)
0094     {
0095       // FIXME
0096       //   - Confirm if TDatabasePDG is needed
0097       //   - Use epsilon comparison here
0098       if (part->get_IonCharge() == 0.0)
0099       {
0100         continue;
0101       }
0102     }
0103
0104     // remove acceptance... _etamin,_etamax
0105     if ((part->get_px() == 0.0) && (part->get_py() == 0.0))
0106     {
0107       continue;  // avoid pt=0
0108     }
0109     float eta = asinh(part->get_pz() / sqrt(pow(part->get_px(), 2) + pow(part->get_py(), 2)));
0110     if (eta < m_EtaMin)
0111     {
0112       continue;
0113     }
0114     if (eta > m_EtaMax)
0115     {
0116       continue;
0117     }
0118
0119     Jet *jet = new Jetv2();
0120     jet->set_px(part->get_px());
0121     jet->set_py(part->get_py());
0122     jet->set_pz(part->get_pz());
0123     jet->set_e(part->get_e());
0124     jet->insert_comp(Jet::PARTICLE, part->get_track_id());
0125     pseudojets.push_back(jet);
0126   }
0127
0128   if (Verbosity() > 0)
0129   {
0130     std::cout << "TruthJetInput::process_event -- exited" << std::endl;
0131   }
0132
0133   return pseudojets;
0134 }