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