g4simulation/g4main/PHG4TruthSubsystem.cc

0001 #include "PHG4TruthSubsystem.h"
0002
0003 #include "PHG4Particle.h"  // for PHG4Particle
0004 #include "PHG4TruthEventAction.h"
0005 #include "PHG4TruthInfoContainer.h"
0006 #include "PHG4TruthTrackingAction.h"
0007
0008 #include <fun4all/Fun4AllReturnCodes.h>
0009
0010 #include <phool/PHCompositeNode.h>
0011 #include <phool/PHIODataNode.h>    // for PHIODataNode
0012 #include <phool/PHNode.h>          // for PHNode
0013 #include <phool/PHNodeIterator.h>  // for PHNodeIterator
0014 #include <phool/PHObject.h>        // for PHObject
0015 #include <phool/getClass.h>
0016 #include <phool/phool.h>  // for PHWHERE
0017
0018 #include <cassert>
0019 #include <cstdlib>  // for exit
0020 #include <iostream>
0021 #include <set>      // for _Rb_tree_iterator, set, _Rb_...
0022 #include <utility>  // for pair
0023
0024 class PHG4EventAction;
0025 class PHG4TrackingAction;
0026
0027 using namespace std;
0028
0029 //_______________________________________________________________________
0030 PHG4TruthSubsystem::PHG4TruthSubsystem(const string& name)
0031   : PHG4Subsystem(name)
0032   , m_EventAction(nullptr)
0033   , m_TrackingAction(nullptr)
0034   , m_SaveOnlyEmbededFlag(false)
0035 {
0036 }
0037
0038 //_______________________________________________________________________
0039 int PHG4TruthSubsystem::InitRun(PHCompositeNode* topNode)
0040 {
0041   PHNodeIterator iter(topNode);
0042   PHCompositeNode* dstNode = dynamic_cast<PHCompositeNode*>(iter.findFirst("PHCompositeNode", "DST"));
0043
0044   // create truth information container
0045   PHG4TruthInfoContainer* truthInfoList = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
0046   if (!truthInfoList)
0047   {
0048     truthInfoList = new PHG4TruthInfoContainer();
0049     dstNode->addNode(new PHIODataNode<PHObject>(truthInfoList, "G4TruthInfo", "PHObject"));
0050   }
0051
0052   // event action
0053   m_EventAction = new PHG4TruthEventAction();
0054
0055   // create tracking action
0056   m_TrackingAction = new PHG4TruthTrackingAction(m_EventAction);
0057
0058   return 0;
0059 }
0060
0061 //_______________________________________________________________________
0062 int PHG4TruthSubsystem::process_event(PHCompositeNode* topNode)
0063 {
0064   // pass top node to event action so that it gets
0065   // relevant nodes needed internally
0066   if (m_EventAction)
0067   {
0068     m_EventAction->SetInterfacePointers(topNode);
0069   }
0070   else
0071   {
0072     cout << PHWHERE << " No EventAction registered" << endl;
0073     exit(1);
0074   }
0075
0076   if (m_TrackingAction)
0077   {
0078     m_TrackingAction->SetInterfacePointers(topNode);
0079   }
0080   else
0081   {
0082     cout << PHWHERE << " No TrackingAction registered" << endl;
0083     exit(1);
0084   }
0085
0086   return Fun4AllReturnCodes::EVENT_OK;
0087 }
0088
0089 int PHG4TruthSubsystem::process_after_geant(PHCompositeNode* topNode)
0090 {
0091   if (m_SaveOnlyEmbededFlag)
0092   {
0093     if (Verbosity() > 1)
0094     {
0095       cout << __PRETTY_FUNCTION__ << " - INFO - only save the G4 truth information that is associated with the embedded particle" << endl;
0096     }
0097
0098     PHG4TruthInfoContainer* truthInfoList = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
0099     assert(truthInfoList);
0100
0101     set<int> savevtxlist;
0102
0103     // remove particle that is not embedd associated
0104     PHG4TruthInfoContainer::Range truth_range = truthInfoList->GetParticleRange();
0105     PHG4TruthInfoContainer::Iterator truthiter = truth_range.first;
0106     while (truthiter != truth_range.second)
0107     {
0108       const int primary_id = (truthiter->second)->get_primary_id();
0109       if (truthInfoList->isEmbeded(primary_id) <= 0)
0110       {
0111         // not a embed associated particle
0112
0113         truthInfoList->delete_particle(truthiter++);
0114       }
0115       else
0116       {
0117         // save vertex id for primary particle which leaves no hit
0118         // in active area
0119         savevtxlist.insert((truthiter->second)->get_vtx_id());
0120         ++truthiter;
0121       }
0122     }
0123
0124     // remove vertex that is not embedd associated
0125     PHG4TruthInfoContainer::VtxRange vtxrange = truthInfoList->GetVtxRange();
0126     PHG4TruthInfoContainer::VtxIterator vtxiter = vtxrange.first;
0127     while (vtxiter != vtxrange.second)
0128     {
0129       if (savevtxlist.find(vtxiter->first) == savevtxlist.end())
0130       {
0131         truthInfoList->delete_vtx(vtxiter++);
0132       }
0133       else
0134       {
0135         ++vtxiter;
0136       }
0137     }
0138   }
0139
0140   return 0;
0141 }
0142
0143 int PHG4TruthSubsystem::ResetEvent(PHCompositeNode* topNode)
0144 {
0145   m_TrackingAction->ResetEvent(topNode);
0146   m_EventAction->ResetEvent(topNode);
0147   return 0;
0148 }
0149
0150 //_______________________________________________________________________
0151 PHG4EventAction* PHG4TruthSubsystem::GetEventAction() const
0152 {
0153   return m_EventAction;
0154 }
0155
0156 //_______________________________________________________________________
0157
0158 PHG4TrackingAction*
0159 PHG4TruthSubsystem::GetTrackingAction() const
0160 {
0161   return m_TrackingAction;
0162 }