g4simulation/g4main/PHG4InEventToTruthInfo.cc

0001 #include "PHG4InEventToTruthInfo.h"
0002
0003 #include "PHG4InEvent.h"
0004 #include "PHG4MCProcessDefs.h"
0005 #include "PHG4Particle.h"
0006 #include "PHG4Particlev2.h"
0007 #include "PHG4Particlev3.h"
0008 #include "PHG4Showerv1.h"
0009 #include "PHG4TruthInfoContainer.h"
0010 #include "PHG4VtxPoint.h"
0011 #include "PHG4VtxPointv2.h"
0012
0013 #include <fun4all/Fun4AllReturnCodes.h>
0014
0015 #include <phool/PHCompositeNode.h>
0016 #include <phool/PHIODataNode.h>
0017 #include <phool/PHNode.h>
0018 #include <phool/PHNodeIterator.h>
0019 #include <phool/PHObject.h>
0020 #include <phool/getClass.h>
0021 #include <phool/phool.h>
0022
0023 #include <TDatabasePDG.h>
0024 #include <TSystem.h>
0025
0026 #include <cmath>
0027 #include <iostream>
0028 #include <iterator>
0029 #include <limits>
0030 #include <map>
0031 #include <memory>
0032
0033 PHG4InEventToTruthInfo::PHG4InEventToTruthInfo(const std::string &name)
0034   : SubsysReco(name)
0035 {
0036 }
0037
0038 int PHG4InEventToTruthInfo::InitRun(PHCompositeNode *topNode)
0039 {
0040   CreateNodeTree(topNode);
0041   return Fun4AllReturnCodes::EVENT_OK;
0042 }
0043
0044 int PHG4InEventToTruthInfo::process_event(PHCompositeNode *topNode)
0045 {
0046   PHG4InEvent *inevent = findNode::getClass<PHG4InEvent>(topNode, "PHG4INEVENT");
0047   if (!inevent)
0048   {
0049     if (Verbosity() > 0)
0050     {
0051       std::cout << PHWHERE << "PHG4INEVENT node not found" << std::endl;
0052     }
0053     return Fun4AllReturnCodes::EVENT_OK;
0054   }
0055
0056   PHG4TruthInfoContainer *truthinfo = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
0057
0058   std::map<int, int> input_to_truth_vtxid;
0059   const auto vtxrange = inevent->GetVertices();
0060   for (auto vtxiter = vtxrange.first; vtxiter != vtxrange.second; ++vtxiter)
0061   {
0062     const int input_vtxid = vtxiter->first;
0063     const int truth_vtxid = truthinfo->maxvtxindex() + 1;
0064     PHG4VtxPoint *truth_vtx = dynamic_cast<PHG4VtxPoint *>(vtxiter->second->CloneMe());
0065     truth_vtx->set_id(truth_vtxid);
0066
0067     const auto inserted = truthinfo->AddVertex(truth_vtxid, truth_vtx);
0068     if (inserted == truthinfo->GetVtxRange().second)
0069     {
0070       std::cout << PHWHERE << " Failure to insert vertex " << truth_vtxid << " into TruthInfo" << std::endl;
0071       truth_vtx->identify();
0072       gSystem->Exit(1);
0073       exit(1);
0074     }
0075
0076     input_to_truth_vtxid[input_vtxid] = truth_vtxid;
0077   }
0078
0079   const auto particlerange = inevent->GetParticles();
0080   // reverse to match how G4 inserts particles into the truth table (last particle first)
0081   for (auto particleiter = std::make_reverse_iterator(particlerange.second);
0082        particleiter != std::make_reverse_iterator(particlerange.first);
0083        ++particleiter)
0084   {
0085     const int input_vtxid = particleiter->first;
0086     PHG4Particle *input_particle = particleiter->second;
0087     if (!input_particle)
0088     {
0089       continue;
0090     }
0091
0092     auto vtxid_iter = input_to_truth_vtxid.find(input_vtxid);
0093     if (vtxid_iter == input_to_truth_vtxid.end())
0094     {
0095       std::cout << PHWHERE << "no vertex " << input_vtxid
0096                 << " found for PHG4InEvent particle" << std::endl;
0097       return Fun4AllReturnCodes::ABORTEVENT;
0098     }
0099
0100     const int truth_trackid = truthinfo->maxtrkindex() + 1;
0101     const int truth_vtxid = vtxid_iter->second;
0102
0103     PHG4Particle *truth_particle = makeParticleCopy(input_particle);
0104     truth_particle->set_track_id(truth_trackid);
0105     truth_particle->set_vtx_id(truth_vtxid);
0106     truth_particle->set_parent_id(0);
0107     truth_particle->set_primary_id(truth_trackid);
0108
0109     const auto inserted = truthinfo->AddParticle(truth_trackid, truth_particle);
0110     if (inserted == truthinfo->GetParticleRange().second)
0111     {
0112       std::cout << PHWHERE << " Failure to insert particle " << truth_trackid << " into TruthInfo" << std::endl;
0113       truth_particle->identify();
0114       gSystem->Exit(1);
0115       exit(1);
0116     }
0117
0118     const int embed_flag = inevent->isEmbeded(input_particle);
0119     if (embed_flag)
0120     {
0121       truthinfo->AddEmbededTrkId(truth_trackid, embed_flag);
0122       truthinfo->AddEmbededVtxId(truth_vtxid, embed_flag);
0123     }
0124   }
0125
0126   return Fun4AllReturnCodes::EVENT_OK;
0127 }
0128
0129 PHG4Particle *PHG4InEventToTruthInfo::makeParticleCopy(PHG4Particle *particle)
0130 {
0131   PHG4Particle *particle_return{nullptr};
0132   if (particle->isIon())
0133   {
0134     particle_return = new PHG4Particlev3(particle);
0135   }
0136   else
0137   {
0138     particle_return = new PHG4Particlev2(particle);
0139   }
0140
0141   if (particle_return->get_name().empty())
0142   {
0143     if (TDatabasePDG::Instance()->GetParticle(particle_return->get_pid()))
0144     {
0145       particle_return->set_name(TDatabasePDG::Instance()->GetParticle(particle_return->get_pid())->GetName());
0146     }
0147     else
0148     {
0149       particle_return->set_name("unknown");
0150     }
0151   }
0152   if (!std::isfinite(particle_return->get_e()))
0153   {
0154     if (TDatabasePDG::Instance()->GetParticle(particle->get_pid()))
0155     {
0156       double m = TDatabasePDG::Instance()->GetParticle(particle->get_pid())->Mass();
0157       double e = sqrt(particle->get_px() * particle->get_px() + particle->get_py() * particle->get_py() + particle->get_pz() * particle->get_pz() + m * m);
0158       particle_return->set_e(e);
0159     }
0160     else
0161     {
0162       particle_return->set_e(std::numeric_limits<double>::quiet_NaN());
0163     }
0164   }
0165   return particle_return;
0166 }
0167
0168 int PHG4InEventToTruthInfo::CreateNodeTree(PHCompositeNode *topNode)
0169 {
0170   PHNodeIterator iter(topNode);
0171   PHCompositeNode *dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
0172   if (!dstNode)
0173   {
0174     std::cout << PHWHERE << "DST node not found" << std::endl;
0175     gSystem->Exit(1);
0176     exit(1);
0177   }
0178   PHG4TruthInfoContainer *truthinfo = findNode::getClass<PHG4TruthInfoContainer>(dstNode, "G4TruthInfo");
0179   if (!truthinfo)
0180   {
0181     truthinfo = new PHG4TruthInfoContainer();
0182     dstNode->addNode(new PHIODataNode<PHObject>(truthinfo, "G4TruthInfo", "PHObject"));
0183   }
0184   return Fun4AllReturnCodes::EVENT_OK;
0185 }