g4simulation/g4jets/TruthJetInput.cc

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