g4simulation/g4main/PHG4TruthTrackingAction.cc

0001 #include "PHG4TruthTrackingAction.h"
0002
0003 #include "PHG4Particle.h"  // for PHG4Particle
0004 #include "PHG4Particlev2.h"
0005 #include "PHG4Particlev3.h"
0006 #include "PHG4ProcessMapPhysics.h"
0007 #include "PHG4Shower.h"  // for PHG4Shower
0008 #include "PHG4Showerv1.h"
0009 #include "PHG4TrackUserInfoV1.h"
0010 #include "PHG4TruthEventAction.h"
0011 #include "PHG4TruthInfoContainer.h"
0012 #include "PHG4UserPrimaryParticleInformation.h"
0013 #include "PHG4VtxPointv2.h"
0014
0015 #include <phool/getClass.h>
0016
0017 #include <Geant4/G4DynamicParticle.hh>     // for G4DynamicParticle
0018 #include <Geant4/G4ParticleDefinition.hh>  // for G4ParticleDefinition
0019 #include <Geant4/G4PrimaryParticle.hh>
0020 #include <Geant4/G4SystemOfUnits.hh>
0021 #include <Geant4/G4ThreeVector.hh>  // for G4ThreeVector
0022 #include <Geant4/G4Track.hh>
0023 #include <Geant4/G4TrackVector.hh>  // for G4TrackVector
0024 #include <Geant4/G4TrackingManager.hh>
0025 #include <Geant4/G4VUserTrackInformation.hh>  // for G4VUserTrackInformation
0026
0027 #include <cmath>     // for sqrt
0028 #include <iostream>  // for operator<<, endl
0029 #include <utility>   // for pair
0030
0031 //________________________________________________________
0032 PHG4TruthTrackingAction::PHG4TruthTrackingAction(PHG4TruthEventAction* eventAction)
0033   : m_EventAction(eventAction)
0034   , m_TruthInfoList(nullptr)
0035   , m_G4ParticleStack()
0036   , m_CurrG4Particle()
0037 {
0038 }
0039
0040 void PHG4TruthTrackingAction::PreUserTrackingAction(const G4Track* track)
0041 {
0042   // insert particle into the output
0043   PHG4Particle* ti = AddParticle(*m_TruthInfoList, *const_cast<G4Track*>(track));
0044
0045   // Negative G4 track id values indicate unwanted tracks to be deleted
0046   // Initially all tracks except primary ones flagged as unwanted
0047   int track_id_g4 = track->GetTrackID() * (track->GetParentID() ? -1 : +1);
0048   int trackid = ti->get_track_id();
0049
0050   // add the user id to the geant4 user info
0051   PHG4TrackUserInfo::SetUserTrackId(const_cast<G4Track*>(track), trackid);
0052
0053   if (PHG4TrackUserInfoV1* p = dynamic_cast<PHG4TrackUserInfoV1*>(track->GetUserInformation()))
0054   {
0055     ti->set_parent_id(p->GetUserParentId());
0056
0057     if (track->GetParentID())
0058     {
0059       ti->set_primary_id(p->GetUserPrimaryId());
0060     }
0061     else
0062     {
0063       PHG4TrackUserInfo::SetUserPrimaryId(const_cast<G4Track*>(track), trackid);
0064       ti->set_primary_id(trackid);
0065     }
0066   }
0067
0068   if (!track->GetParentID())
0069   {
0070     // primary track - propagate the barcode information
0071     PHG4UserPrimaryParticleInformation* userdata = static_cast<PHG4UserPrimaryParticleInformation*>(track->GetDynamicParticle()->GetPrimaryParticle()->GetUserInformation());
0072     if (userdata)
0073     {
0074       ti->set_barcode(userdata->get_user_barcode());
0075     }
0076   }
0077
0078   int vtxindex = ti->get_vtx_id();
0079
0080
0081   // maybe we should do the sPHENIX primary tracking here as here is the place where the parent id etc. are finally set
0082
0083   if (issPHENIXPrimary(*m_TruthInfoList, ti))
0084   {
0085
0086     //we also want to set keep this track
0087     PHG4TrackUserInfoV1* userinfo = dynamic_cast<PHG4TrackUserInfoV1*>(track->GetUserInformation());
0088     if (userinfo)
0089     {
0090       userinfo->SetKeep(true);
0091     }
0092
0093     PHG4Particle* newparticle = dynamic_cast<PHG4Particle*>(ti->CloneMe());
0094
0095     m_TruthInfoList->AddsPHENIXPrimaryParticle(trackid, newparticle);
0096   }
0097
0098   m_CurrG4Particle = {track_id_g4, trackid, vtxindex};
0099
0100   // create or add to a new shower object --------------------------------------
0101   if (!track->GetParentID())
0102   {
0103     PHG4Showerv1* shower = new PHG4Showerv1();
0104     PHG4TrackUserInfo::SetShower(const_cast<G4Track*>(track), shower);
0105     m_TruthInfoList->AddShower(trackid, shower);
0106     shower->set_id(trackid);  // fyi, secondary showers may not share these ids
0107     shower->set_parent_particle_id(trackid);
0108     shower->set_parent_shower_id(0);
0109   }
0110   else
0111   {
0112     // get shower
0113     if (G4VUserTrackInformation* p = track->GetUserInformation())
0114     {
0115       if (PHG4TrackUserInfoV1* pp = dynamic_cast<PHG4TrackUserInfoV1*>(p))
0116       {
0117         if (pp->GetShower())
0118         {
0119           pp->GetShower()->add_g4particle_id(trackid);
0120           pp->GetShower()->add_g4vertex_id(vtxindex);
0121         }
0122       }
0123     }
0124   }
0125
0126   // tell the primary particle copy in G4 where this output will be stored
0127   if (!track->GetParentID())
0128   {
0129     PHG4UserPrimaryParticleInformation* userdata = static_cast<PHG4UserPrimaryParticleInformation*>(track->GetDynamicParticle()->GetPrimaryParticle()->GetUserInformation());
0130     if (userdata)
0131     {
0132       userdata->set_user_track_id(trackid);
0133       userdata->set_user_vtx_id(vtxindex);
0134     }
0135   }
0136
0137   return;
0138 }
0139
0140 void PHG4TruthTrackingAction::PostUserTrackingAction(const G4Track* track)
0141 {
0142   if (fpTrackingManager)
0143   {
0144     int trackid = track->GetTrackID();
0145     int primaryid = 0;
0146     PHG4Shower* shower = nullptr;
0147     if (PHG4TrackUserInfoV1* p = dynamic_cast<PHG4TrackUserInfoV1*>(track->GetUserInformation()))
0148     {
0149       trackid = p->GetUserTrackId();
0150       primaryid = p->GetUserPrimaryId();
0151       shower = p->GetShower();
0152     }
0153
0154     G4TrackVector* secondaries = fpTrackingManager->GimmeSecondaries();
0155     if (secondaries)
0156     {
0157       for (auto secondary : *secondaries)
0158       {
0159         PHG4TrackUserInfo::SetUserParentId(const_cast<G4Track*>(secondary), trackid);
0160         PHG4TrackUserInfo::SetUserPrimaryId(const_cast<G4Track*>(secondary), primaryid);
0161         PHG4TrackUserInfo::SetShower(const_cast<G4Track*>(secondary), shower);
0162       }
0163     }
0164   }
0165
0166   if (PHG4TrackUserInfoV1* p = dynamic_cast<PHG4TrackUserInfoV1*>(track->GetUserInformation()))
0167   {
0168     if (p->GetKeep())
0169     {
0170       int trackid = p->GetUserTrackId();
0171       m_EventAction->AddTrackidToWritelist(trackid);
0172     }
0173   }
0174
0175   UpdateG4ParticleStack(track);
0176 }
0177
0178 /**
0179  * Updates the stack of parent particles and removes unwanted ones from the
0180  * truth info container.
0181  */
0182 void PHG4TruthTrackingAction::UpdateG4ParticleStack(const G4Track* track)
0183 {
0184   while (!m_G4ParticleStack.empty())
0185   {
0186     if (std::abs(m_G4ParticleStack.back().g4track_id) == track->GetParentID())
0187     {
0188       break;
0189     }
0190     else
0191     {
0192       if (m_G4ParticleStack.back().g4track_id < 0)
0193       {
0194         m_TruthInfoList->delete_particle(m_G4ParticleStack.back().particle_id);
0195       }
0196       m_G4ParticleStack.pop_back();
0197     }
0198   }
0199
0200   m_G4ParticleStack.push_back(m_CurrG4Particle);
0201
0202   // Change sign of G4 track id of all upstream tracks in the stack to positive
0203   // in order to keep the track
0204   PHG4TrackUserInfoV1* p = dynamic_cast<PHG4TrackUserInfoV1*>(track->GetUserInformation());
0205   bool keep_curr_track = p && p->GetKeep() ? true : false;
0206
0207   auto stack_iter = m_G4ParticleStack.rbegin();
0208   while (keep_curr_track && stack_iter != m_G4ParticleStack.rend() && stack_iter->g4track_id < 0)
0209   {
0210     stack_iter->g4track_id *= -1;
0211     ++stack_iter;
0212   }
0213 }
0214
0215 void PHG4TruthTrackingAction::SetInterfacePointers(PHCompositeNode* topNode)
0216 {
0217   // now look for the map and grab a pointer to it.
0218   m_TruthInfoList = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
0219
0220   // if we do not find the node we need to make it.
0221   if (!m_TruthInfoList)
0222   {
0223     std::cout << "PHG4TruthEventAction::SetInterfacePointers - unable to find G4TruthInfo" << std::endl;
0224   }
0225 }
0226
0227 int PHG4TruthTrackingAction::ResetEvent(PHCompositeNode* /*unused*/)
0228 {
0229   m_VertexMap.clear();
0230
0231   while (!m_G4ParticleStack.empty())
0232   {
0233     if (m_G4ParticleStack.back().g4track_id < 0)
0234     {
0235       m_TruthInfoList->delete_particle(m_G4ParticleStack.back().particle_id);
0236     }
0237     m_G4ParticleStack.pop_back();
0238   }
0239
0240   return 0;
0241 }
0242
0243 PHG4Particle* PHG4TruthTrackingAction::AddParticle(PHG4TruthInfoContainer& truth, G4Track& track)
0244 {
0245   int trackid = 0;
0246   if (track.GetParentID())
0247   {
0248     // secondaries get negative user ids and increment downward between geant subevents
0249     trackid = truth.mintrkindex() - 1;
0250   }
0251   else
0252   {
0253     // primaries get positive user ids and increment upward between geant subevents
0254     trackid = truth.maxtrkindex() + 1;
0255   }
0256
0257   // determine the momentum vector
0258   G4ParticleDefinition* def = track.GetDefinition();
0259   int pdgid = def->GetPDGEncoding();
0260   double mass = def->GetPDGMass();
0261   double ke = track.GetVertexKineticEnergy();
0262   double ptot = sqrt(ke * ke + 2.0 * mass * ke);
0263   G4ThreeVector pdir = track.GetVertexMomentumDirection();
0264   pdir *= ptot;
0265   PHG4Particle* ti = nullptr;
0266   // create a new particle -----------------------------------------------------
0267   if (def->IsGeneralIon())  // for ions save a and z in v3 of phg4particle
0268   {
0269     ti = new PHG4Particlev3();
0270     ti->set_A(def->GetAtomicMass());
0271     ti->set_Z(def->GetAtomicNumber());
0272   }
0273   else
0274   {
0275     ti = new PHG4Particlev2;
0276   }
0277   ti->set_px(pdir[0] / GeV);
0278   ti->set_py(pdir[1] / GeV);
0279   ti->set_pz(pdir[2] / GeV);
0280   ti->set_track_id(trackid);
0281
0282   ti->set_parent_id(track.GetParentID());
0283   ti->set_primary_id(trackid);
0284
0285   ti->set_pid(pdgid);
0286   ti->set_name(def->GetParticleName());
0287   ti->set_e(track.GetTotalEnergy() / GeV);
0288
0289   // Add new or reuse a vertex and let the track know about it
0290   PHG4VtxPoint* vtx = AddVertex(truth, track);
0291   ti->set_vtx_id(vtx->get_id());
0292
0293   // use a new map to hold the new primary particle list
0294   // some debug print
0295   /*
0296   {
0297     float vtx_x = vtx->get_x() * cm;
0298     float vtx_y = vtx->get_y() * cm;
0299     float vtx_z = vtx->get_z() * cm;
0300     auto process = vtx->get_process();
0301     std::cout << "PHG4TruthTrackingAction::AddParticle - Adding particle with track id " << trackid
0302               << ", vtx id " << ti->get_vtx_id()
0303               << ", vtx position (" << vtx_x << ", " << vtx_y << ", " << vtx_z << ") cm"
0304               << ", process: " << process
0305               << ", parent id: " << ti->get_parent_id()
0306               << ", pid: " << ti->get_pid()
0307               << ", name: " << ti->get_name()
0308               << std::endl;
0309   }
0310   */
0311   return truth.AddParticle(trackid, ti)->second;
0312 }
0313
0314 PHG4VtxPoint* PHG4TruthTrackingAction::AddVertex(PHG4TruthInfoContainer& truth, const G4Track& track)
0315 {
0316   G4ThreeVector v = track.GetVertexPosition();
0317   int vtxindex = (track.GetParentID() == 0 ? truth.maxvtxindex() + 1 : truth.minvtxindex() - 1);
0318
0319   auto [iter, inserted] = m_VertexMap.insert(std::make_pair(v, vtxindex));
0320
0321   // If could not add a unique vertex => return the existing one
0322   if (!inserted)
0323   {
0324     return truth.GetVtxMap().find(iter->second)->second;
0325   }
0326   // get G4Track creator process
0327   const auto g4Process = track.GetCreatorProcess();
0328   // convert G4 Process to MC process
0329   const auto process = PHG4ProcessMapPhysics::Instance().GetMCProcess(g4Process);
0330   // otherwise, create and add a new one
0331   PHG4VtxPoint* vtxpt = new PHG4VtxPointv2(v[0] / cm, v[1] / cm, v[2] / cm, track.GetGlobalTime() / ns, vtxindex, process);
0332
0333   return truth.AddVertex(vtxindex, vtxpt)->second;
0334 }
0335
0336 bool PHG4TruthTrackingAction::issPHENIXPrimary(PHG4TruthInfoContainer& truth, PHG4Particle* particle) const
0337 {
0338   PHG4VtxPoint* vtx = truth.GetVtx(particle->get_vtx_id());
0339   if (!vtx)
0340   {
0341     // something is very very wrong... I guess
0342     std::cerr << "PHG4TruthTrackingAction::issPHENIXPrimary - no vertex found for particle with track id " << particle->get_track_id() << std::endl;
0343     return false;
0344   }
0345   auto process = vtx->get_process();
0346   int pdgid = particle->get_pid();
0347   // if not long-lived, then it is not a primary
0348   if (!isLongLived(pdgid))
0349   {
0350     return false;
0351   }
0352   // check the production process
0353   // if not decay or primary, then it is not a primary
0354   // debug print for pid, track id, parent id, and process
0355   /*
0356   std::cout << "PHG4TruthTrackingAction::issPHENIXPrimary - checking particle with track id " << particle->get_track_id()
0357             << ", pid: " << pdgid
0358             << ", parent id: " << particle->get_parent_id()
0359             << ", process: " << process
0360             << std::endl;
0361             */
0362   if (!(process == PHG4MCProcess::kPPrimary || process == PHG4MCProcess::kPDecay) && particle->get_parent_id())  // all primary particles seems to have unkown process id
0363   {
0364     return false;
0365   }
0366   // now we are clear from particle produced from material interactions
0367   if (particle->get_parent_id() == 0)
0368   {
0369     // conditioning on the above, if the track is primary, then it is a sPHENIX primary
0370     return true;
0371   }
0372   // not we want to check if their parent is long-lived or primary
0373   // in G4 parent should always process before child, so we can just go up the tree
0374   PHG4Particle* parent = truth.GetParticle(particle->get_parent_id());
0375   // if there is a loop of the parent, then we have a problem lol with this while loop btw
0376   while (parent)
0377   {
0378     if (isLongLived(parent->get_pid()))
0379     {
0380       // if the parent is long-lived or primary, then it is not a sPHENIX primary
0381       return false;
0382     }
0383     PHG4VtxPoint* vtx_parent = truth.GetVtx(parent->get_vtx_id());
0384     if (!vtx_parent)
0385     {
0386       // something is very very wrong... I guess
0387       std::cerr << "PHG4TruthTrackingAction::issPHENIXPrimary - no vertex found for parent particle with track id " << parent->get_track_id() << std::endl;
0388       return false;
0389     }
0390     process = vtx_parent->get_process();
0391     // if parent is not from decay or primary, then it is not a sPHENIX primary
0392     if (!(process == PHG4MCProcess::kPPrimary || process == PHG4MCProcess::kPDecay) && parent->get_parent_id())
0393     {
0394       return false;
0395     }
0396     // otherwise, go up the tree
0397     parent = truth.GetParticle(parent->get_parent_id());
0398   }
0399
0400   return true;
0401 }
0402
0403 bool PHG4TruthTrackingAction::isLongLived(int pid) const
0404 {
0405   // see https://inspirehep.net/files/4c26ef5fb432df99bdc1ff847653502f
0406   // Check nuclus
0407   if (pid > 1000000000) return true;
0408   // this needs to be hardcoded somehow... :(
0409   // but in the future we can find a better home for this piece of code
0410   switch (pid)
0411   {
0412   case 11:     // electron
0413   case -11:    // positron
0414   case 13:     // muon
0415   case -13:    // antimuon
0416   case 22:     // photon
0417   case 211:    // pi+
0418   case -211:   // pi-
0419   case 321:    // K+
0420   case -321:   // K-
0421   case 310:    // K0S
0422   case 130:    // K0L
0423   case 2212:   // proton
0424   case -2212:  // antiproton
0425   case 2112:   // neutron
0426   case -2112:  // antineutron
0427   case 3122:   // Lambda
0428   case -3122:  // anti-Lambda
0429   case 3222:   // Sigma+
0430   case -3222:  // anti-Sigma+ ?
0431   case 3112:   // Sigma-
0432   case -3112:  // anti-Sigma-
0433   case 3312:   // Xi-
0434   case -3312:  // anti-Xi-
0435   case 3322:   // Xi0
0436   case -3322:  // anti-Xi0
0437   case 3334:   // Omega-
0438   case -3334:  // anti-Omega-
0439   case 12:     // neutrino
0440   case -12:    // antineutrino
0441   case 14:     // muon neutrino
0442   case -14:    // muon antineutrino
0443   case 16:     // tau neutrino
0444   case -16:    // tau antineutrino
0445
0446     return true;
0447   default:
0448     return false;
0449   }
0450 }