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