tca::MCParticleListUtils Class Reference

#include "TCTruth.h"

Public Member Functions
	MCParticleListUtils (TjStuff &my_tjs)
void	MakeTruTrajPoint (unsigned int MCParticleListIndex, TrajPoint &tp)
ShowerStruct3D	MakeCheatShower (unsigned int mcpIndex, Point3_t primVx, int &truParentPFP)
bool	PrimaryElectronStart (Point3_t &start, Vector3_t &dir, float &energy)
int	PrimaryElectronPFPID (const geo::TPCID &tpcid)
int	PrimaryElectronTjID (CTP_t inCTP)
int	MCParticleStartTjID (unsigned int MCParticleListIndex, CTP_t inCTP)
unsigned int	GetMCPartListIndex (const TrajPoint &tp)
unsigned int	GetMCPartListIndex (const Trajectory &tj, unsigned short &nTruHits)
unsigned int	GetMCPartListIndex (const ShowerStruct &ss, unsigned short &nTruHits)

Public Attributes
TjStuff &	tjs

Detailed Description

Definition at line 81 of file TCTruth.h.

Constructor & Destructor Documentation

tca::MCParticleListUtils::MCParticleListUtils ( TjStuff &  my_tjs )

inline

Definition at line 85 of file TCTruth.h.

: tjs(my_tjs) {}

Member Function Documentation

unsigned int tca::MCParticleListUtils::GetMCPartListIndex

(

const TrajPoint &

tp

)

Definition at line 1499 of file TCTruth.cxx.

References tca::TrajPoint::Chg, tca::TjStuff::fHits, tca::TrajPoint::Hits, tca::TjStuff::MCPartList, tca::TruthMatcher::tjs, and tca::TrajPoint::UseHit.

  {
    // Returns the MCParticle index that best matches the hits used in the Tp
    if(tjs.MCPartList.empty()) return UINT_MAX;
    if(tp.Chg <= 0) return UINT_MAX;
    std::vector<unsigned int> mcpIndex;
    std::vector<unsigned short> mcpCnt;
    for(unsigned short ii = 0; ii < tp.Hits.size(); ++ii) {
      if(!tp.UseHit[ii]) continue;
      unsigned int mcpi = tjs.fHits[tp.Hits[ii]].MCPartListIndex;
      if(mcpi == UINT_MAX) continue;
      unsigned short indx = 0;
      for(indx = 0; indx < mcpIndex.size(); ++indx) if(mcpi == mcpIndex[indx]) break;
      if(indx == mcpIndex.size()) {
        // not in the list so add it
        mcpIndex.push_back(mcpi);
        mcpCnt.push_back(1);
      } else {
        ++mcpCnt[indx];
      }
    } // ii
    if(mcpIndex.empty()) return UINT_MAX;
    if(mcpIndex.size() == 1) return mcpIndex[0];
    unsigned int indx = 0;
    unsigned short maxCnt = 0;
    for(unsigned short ii = 0; ii < mcpIndex.size(); ++ii) {
      if(mcpCnt[ii] > maxCnt) {
        maxCnt = mcpCnt[ii];
        indx = mcpIndex[ii];
      }
    } // ii
    return indx;
  } // GetMCPartListIndex

unsigned int tca::MCParticleListUtils::GetMCPartListIndex	(	const Trajectory &	tj,
		unsigned short &	nTruHits
	)

Definition at line 1571 of file TCTruth.cxx.

References tca::TjStuff::fHits, tca::TjStuff::MCPartList, tca::Trajectory::Pts, and tca::TruthMatcher::tjs.

  {
    // Returns the index of the MCParticle that has the most number of matches
    // to the hits in this trajectory
    
    if(tjs.MCPartList.empty()) return UINT_MAX;
    
    // Check all hits associated with this Tj
    std::vector<unsigned int> pListCnt(tjs.MCPartList.size());
    
    for(auto& tp : tj.Pts) {
      for(unsigned short ii = 0; ii < tp.Hits.size(); ++ii) {
        if(!tp.UseHit[ii]) continue;
        unsigned int iht = tp.Hits[ii];
        // ignore unmatched hits
        if(tjs.fHits[iht].MCPartListIndex > tjs.MCPartList.size() - 1) continue;
        ++pListCnt[tjs.fHits[iht].MCPartListIndex];
      } // ii
    } // pt
    
    unsigned int pIndex = UINT_MAX;
    nTruHits = 0;
    for(unsigned short ii = 0; ii < pListCnt.size(); ++ii) {
      if(pListCnt[ii] > nTruHits) {
        nTruHits = pListCnt[ii];
        pIndex = ii;
      }
    } // ii
    
    return pIndex;
    
  } // GetMCPartListIndex

unsigned int tca::MCParticleListUtils::GetMCPartListIndex	(	const ShowerStruct &	ss,
		unsigned short &	nTruHits
	)

Definition at line 1534 of file TCTruth.cxx.

References tca::TjStuff::allTraj, tca::TjStuff::fHits, tca::TjStuff::MCPartList, tca::Trajectory::Pts, tca::ShowerStruct::TjIDs, and tca::TruthMatcher::tjs.

  {
    // Returns the index of the MCParticle that has the most number of matches
    // to the hits in this shower
    
    if(tjs.MCPartList.empty()) return UINT_MAX;
    if(ss.TjIDs.empty()) return UINT_MAX;
    
    std::vector<unsigned int> pListCnt(tjs.MCPartList.size());
    
    for(auto& tjid : ss.TjIDs) {
      Trajectory& tj = tjs.allTraj[tjid - 1];
      for(auto& tp : tj.Pts) {
        for(unsigned short ii = 0; ii < tp.Hits.size(); ++ii) {
          if(!tp.UseHit[ii]) continue;
          unsigned int iht = tp.Hits[ii];
          // ignore unmatched hits
          if(tjs.fHits[iht].MCPartListIndex > tjs.MCPartList.size() - 1) continue;
          ++pListCnt[tjs.fHits[iht].MCPartListIndex];
        } // ii
      } // pt
    } // tjid
    
    unsigned int pIndex = UINT_MAX;
    nTruHits = 0;
    for(unsigned short ii = 0; ii < pListCnt.size(); ++ii) {
      if(pListCnt[ii] > nTruHits) {
        nTruHits = pListCnt[ii];
        pIndex = ii;
      }
    } // ii
    
    return pIndex;
    
  } // GetMCPartListIndex

ShowerStruct3D tca::MCParticleListUtils::MakeCheatShower	(	unsigned int	mcpIndex,
		Point3_t	primVx,
		int &	truParentPFP
	)

Definition at line 1191 of file TCTruth.cxx.

References tca::TjStuff::allTraj, close(), tca::CreateSS(), tca::CreateSS3(), geo::CryostatID::Cryostat, tca::EncodeCTP(), evd::details::end(), sim::ParticleList::EveId(), tca::TjStuff::fHits, tca::GetAssns(), tca::kMat3D, tca::TjStuff::MCPartList, tca::TjStuff::NumPlanes, cheat::ParticleInventoryService::ParticleList(), tca::TjStuff::pfps, tca::PosSep(), ss, tca::StoreShower(), tca::TruthMatcher::tjs, geo::TPCID::TPC, tca::TjStuff::TPCID, tca::TjStuff::UnitsPerTick, tca::UpdateShower(), and tca::TjStuff::WirePitch.

  {
    // Make a 3D shower for an electron or photon MCParticle in the current TPC. 
    // The shower ID is set to 0 if there is a failure. The ID of the most likely parent pfp is returned
    // as well - the one that is closest to the primary vertex position
    auto ss3 = CreateSS3(tjs, tjs.TPCID);
    truParentPFP = 0;
    // save the ID 
    int goodID = ss3.ID;
    // set it to the failure state
    ss3.ID = 0;
    ss3.Cheat = true;
    if(mcpIndex > tjs.MCPartList.size() - 1) return ss3;
    auto& mcp = tjs.MCPartList[mcpIndex];
    int pdg = abs(mcp->PdgCode());
    if(!(pdg == 11 || pdg == 111)) return ss3;
    art::ServiceHandle<cheat::ParticleInventoryService> pi_serv;
    int eveID = mcp->TrackId();
    float showerEnergy = 1000 * mcp->E();
    float shMaxAlong = 7.0 * log(showerEnergy / 15);
    std::cout<<"MCS: showerEnergy "<<(int)showerEnergy<<" MeV. shMaxAlong "<<shMaxAlong<<" cm\n";
    
    // put the shower hits in two vectors. One for the InShower hits
    // and one for the shower parent
    std::vector<std::vector<unsigned int>> showerHits(tjs.NumPlanes);
    std::vector<std::vector<unsigned int>> parentHits(tjs.NumPlanes);
    unsigned short nsh = 0;
    unsigned short npar = 0;
    unsigned int cstat = tjs.TPCID.Cryostat;
    unsigned int tpc = tjs.TPCID.TPC;
    for(unsigned int iht = 0; iht < tjs.fHits.size(); ++iht) {
      auto& hit = tjs.fHits[iht];
      if(hit.MCPartListIndex > tjs.MCPartList.size()-1) continue;
      if(hit.ArtPtr->WireID().TPC != tpc) continue;
      if(hit.ArtPtr->WireID().Cryostat != cstat) continue;
      if(hit.Integral <= 0) continue;
      // require that it be used in a tj to ignore hits that are far
      // from the shower
      if(hit.InTraj <= 0) continue;
      auto& shmcp = tjs.MCPartList[hit.MCPartListIndex];
      // look for an electron
      if(abs(shmcp->PdgCode()) != 11) continue;
      // with eve ID = the primary electron being considered
      int eid = pi_serv->ParticleList().EveId(shmcp->TrackId());
      if(eid != eveID) continue;
      if(shmcp->TrackId() == eid) {
        // store the shower parent hit
        parentHits[hit.ArtPtr->WireID().Plane].push_back(iht);
        ++npar;
      } else {
        // store the InShower hit
        showerHits[hit.ArtPtr->WireID().Plane].push_back(iht);
        ++nsh;
      }
    } // iht
    // require more hits in the shower than in the parent
    if(npar > nsh) return ss3;
/*
    for(unsigned short plane = 0; plane < tjs.NumPlanes; ++plane) {
      std::cout<<" plane "<<plane<<" shower hits "<<showerHits[plane].size();
      std::cout<<" parentHits "<<parentHits[plane].size()<<"\n";
    } // plane
*/
    // create 2D cheat showers in each plane
    std::vector<int> dummy_tjlist;
    for(unsigned short plane = 0; plane < tjs.NumPlanes; ++plane) {
      CTP_t inCTP = EncodeCTP(cstat, tpc, plane);
      auto ss = CreateSS(tjs, inCTP, dummy_tjlist);
      // fill the ShPts
      ss.ShPts.resize(showerHits[plane].size());
      for(unsigned short iht = 0; iht < showerHits[plane].size(); ++iht) {
        auto& hit = tjs.fHits[showerHits[plane][iht]];
        ss.ShPts[iht].HitIndex = showerHits[plane][iht];
        ss.ShPts[iht].TID = 0;
        ss.ShPts[iht].Chg = hit.Integral;
        ss.ShPts[iht].Pos[0] = hit.ArtPtr->WireID().Wire;
        ss.ShPts[iht].Pos[1] = hit.PeakTime * tjs.UnitsPerTick;
      } // iht
      ss.SS3ID = goodID;
      if(!UpdateShower("MCS", tjs, ss, true)) {
        std::cout<<"Failed to update 2S"<<ss.ID<<"\n";
        return ss3;
      } // UpdateShower failed
      // remove shower points that are far away
      std::vector<ShowerPoint> spts;
      for(auto& spt : ss.ShPts) {
        // don't make a tight cut right now. The shower parameterization isn't great
        double along = tjs.WirePitch * spt.RotPos[0];
        float tau = along / shMaxAlong;
//        auto& hit = tjs.fHits[spt.HitIndex];
//        if(ss.ID == 2) std::cout<<"chk "<<PrintHit(hit)<<" along "<<(int)along<<" tau "<<std::fixed<<std::setprecision(1)<<tau<<"\n";
        if(tau > -1 && tau < 2) {
          spts.push_back(spt);
        } else {
//          std::cout<<"Skip "<<PrintHit(hit)<<" along "<<(int)along<<"\n";
        }
      } // spt
      std::cout<<" 2S"<<ss.ID<<" shpts "<<ss.ShPts.size()<<" new "<<spts.size()<<"\n";
      ss.ShPts = spts;
      ss.NeedsUpdate = true;
      if(!UpdateShower("MCS", tjs, ss, true)) {
        std::cout<<"Failed to update 2S"<<ss.ID<<"\n";
        return ss3;
      } // UpdateShower failed
      if(!StoreShower("MCS", tjs, ss)) {
        std::cout<<"Failed to store 2S"<<ss.ID<<"\n";
        return ss3;
      } // UpdateShower failed
      ss3.CotIDs.push_back(ss.ID);
    } // plane
    ss3.ID = goodID;
    if(!UpdateShower("MCS", tjs, ss3, true)) {
      std::cout<<"SS3 Failed...\n";
      ss3.ID = 0;
      return ss3;
    }
    // success
    
    // now look for a parent pfp
    std::vector<std::pair<int, int>> pcnt;
    for(unsigned short plane = 0; plane < tjs.NumPlanes; ++plane) {
      for(auto iht : parentHits[plane]) {
        auto& hit = tjs.fHits[iht];
        if(hit.InTraj <= 0) continue;
        auto& tj = tjs.allTraj[hit.InTraj - 1];
        if(!tj.AlgMod[kMat3D]) continue;
        auto TInP = GetAssns(tjs, "T", tj.ID, "P");
        if(TInP.empty()) continue;
        auto& pfp = tjs.pfps[TInP[0] - 1];
        unsigned short indx = 0;
        for(indx = 0; indx < pcnt.size(); ++indx) if(pfp.ID == pcnt[indx].first) break;
        if(indx == pcnt.size()) {
          pcnt.push_back(std::make_pair(pfp.ID, 1));
        } else {
          ++pcnt[indx].second;
        }
      } // iht
    } // plane
    float close = 5;
    for(auto pcn : pcnt) {
      if(pcn.second < 6) continue;
      auto& pfp = tjs.pfps[pcn.first - 1];
      for(unsigned short end = 0; end < 2; ++end) {
        float sep = PosSep(pfp.XYZ[end], primVx);
        if(sep > close) continue;
        close = sep;
        truParentPFP = pfp.ID;
      }
    } // pcn
    std::cout<<"truParent P"<<truParentPFP<<" close "<<std::fixed<<std::setprecision(2)<<close<<"\n";
    
    return ss3;
  } // MakeCheatShower

void tca::MCParticleListUtils::MakeTruTrajPoint	(	unsigned int	MCParticleListIndex,
		TrajPoint &	tp
	)

Definition at line 1174 of file TCTruth.cxx.

References tca::TrajPoint::CTP, dir, tca::MakeBareTP(), tca::TjStuff::MCPartList, simb::MCParticle::Px(), simb::MCParticle::Py(), simb::MCParticle::Pz(), tca::SetMag(), tca::TruthMatcher::tjs, simb::MCParticle::Vx(), simb::MCParticle::Vy(), and simb::MCParticle::Vz().

  {
    // Creates a trajectory point at the start of the MCParticle with index MCParticleListIndex. The
    // projected length of the MCParticle in the plane coordinate system is stored in TruTp.Delta.
    // The calling function should specify the CTP in which this TP resides.
    
    if(MCParticleListIndex > tjs.MCPartList.size() - 1) return;
    
    const simb::MCParticle* mcp = tjs.MCPartList[MCParticleListIndex];
    
    Point3_t pos {{mcp->Vx(), mcp->Vy(), mcp->Vz()}};
    Vector3_t dir {{mcp->Px(), mcp->Py(), mcp->Pz()}};
    SetMag(dir, 1);
    tp = MakeBareTP(tjs, pos, dir, tp.CTP);
  } // MakeTruTrajPoint

int tca::MCParticleListUtils::MCParticleStartTjID	(	unsigned int	MCParticleListIndex,
		CTP_t	inCTP
	)

Definition at line 1454 of file TCTruth.cxx.

References geo::CryostatID::Cryostat, tca::DecodeCTP(), tca::TjStuff::fHits, tca::TjStuff::MCPartList, geo::PlaneID::Plane, tca::TruthMatcher::tjs, tmp, and geo::TPCID::TPC.

  {
    // Finds the trajectory that has hits matched to the MC Particle and is the closest to the
    // MCParticle start vertex
    
    if(mcpIndex > tjs.MCPartList.size() - 1) return 0;
    
    geo::PlaneID planeID = DecodeCTP(inCTP);
    
    // tj ID and occurrence count
    std::vector<std::array<int, 2>> t_cnt;
    for(auto hit : tjs.fHits) {
      if(hit.InTraj <= 0) continue;
      if(hit.MCPartListIndex != mcpIndex) continue;
      if(hit.ArtPtr->WireID().TPC != planeID.TPC) continue;
      if(hit.ArtPtr->WireID().Plane != planeID.Plane) continue;
      if(hit.ArtPtr->WireID().Cryostat != planeID.Cryostat) continue;
      unsigned short indx = 0;
      for(indx = 0; indx < t_cnt.size(); ++indx) if(t_cnt[indx][0] == hit.InTraj) break;
      if(indx == t_cnt.size()) {
        // didn't find the traj ID in t_cnt so add it
        std::array<int, 2> tmp;
        tmp[0] = hit.InTraj;
        tmp[1] = 1;
        t_cnt.push_back(tmp);
      } else {
        // count occurrences of this tj ID
        ++t_cnt[indx][1];
      }
    } // hit
    
    if(t_cnt.empty()) return 0;
    unsigned short occMax = 0;
    unsigned short occIndx = 0;
    for(unsigned short ii = 0; ii < t_cnt.size(); ++ii) {
      auto& tcnt = t_cnt[ii];
      if(tcnt[1] < occMax) continue;
      occMax = tcnt[1];
      occIndx = ii;
    } // tcnt
    return t_cnt[occIndx][0];
    
  } // MCParticleStartTj

int tca::MCParticleListUtils::PrimaryElectronPFPID

(

const geo::TPCID &

tpcid

)

Definition at line 1374 of file TCTruth.cxx.

References tca::TjStuff::allTraj, close(), geo::CryostatID::Cryostat, tca::EncodeCTP(), sim::ParticleList::EveId(), tca::FarEnd(), tca::TjStuff::fHits, tca::GetAssns(), tca::kMat3D, tca::kUsedHits, tca::TjStuff::MCPartList, tca::TjStuff::NumPlanes, part, cheat::ParticleInventoryService::ParticleList(), tca::TjStuff::pfps, tca::PosSep2(), tca::PutTrajHitsInVector(), tca::SetIntersection(), tca::TruthMatcher::tjs, and geo::TPCID::TPC.

  {
    // Returns the ID of a pfp that has hits whose eve ID is a primary electron
    // and is the closest (< 5 WSE units) to the primary electron start
    if(tjs.MCPartList.empty()) return 0;
    art::ServiceHandle<cheat::ParticleInventoryService> pi_serv;
    TruthMatcher tm{tjs};
    for(unsigned int part = 0; part < tjs.MCPartList.size(); ++part) {
      auto& mcp = tjs.MCPartList[part];
      // require electron
      if(abs(mcp->PdgCode()) != 11) continue;
      int eveID = pi_serv->ParticleList().EveId(mcp->TrackId());
      if(mcp->TrackId() != eveID) continue;
      Point3_t start = {{mcp->Vx(), mcp->Vy(), mcp->Vz()}};
      auto const* SCE = lar::providerFrom<spacecharge::SpaceChargeService>();
      geo::Vector_t posOffsets = SCE->GetPosOffsets({start[0], start[1], start[2]});
      posOffsets.SetX(-posOffsets.X());
      start[0] += posOffsets.X();
      start[1] += posOffsets.Y();
      start[2] += posOffsets.Z();
      // make a list of pfps that use tjs that use these hits
      std::vector<int> pfplist;
      for(unsigned short plane = 0; plane < tjs.NumPlanes; ++plane) {
        CTP_t inCTP = EncodeCTP(tpcid.Cryostat, tpcid.TPC, plane);
        std::vector<unsigned int> mcphits = tm.PutMCPHitsInVector(part, inCTP);
        if(mcphits.empty()) continue;
        for(auto iht : mcphits) {
          auto& hit = tjs.fHits[iht];
          if(hit.InTraj <= 0) continue;
          // require that the tj is 3D-matched
          auto& tj = tjs.allTraj[hit.InTraj - 1];
          if(!tj.AlgMod[kMat3D]) continue;
          // find the number of hits that are in mcphits
          auto tjhits = PutTrajHitsInVector(tj, kUsedHits);
          auto shared = SetIntersection(tjhits, mcphits);
          if(shared.size() < 10) continue;
          // find out what pfp it is used in.
          auto TInP = GetAssns(tjs, "T", tj.ID, "P");
          if(TInP.size() != 1) continue;
          int pid = TInP[0];
          if(std::find(pfplist.begin(), pfplist.end(), pid) == pfplist.end()) pfplist.push_back(pid);
        } // iht
      } // plane
      if(pfplist.empty()) return 0;
      // Use the one that is closest to the true start position, not the
      // one that has the most matching hits. Electrons are likely to be
      // poorly reconstructed
      int pfpid = 0;
      float close = 1E6;
      for(auto pid : pfplist) {
        auto& pfp = tjs.pfps[pid - 1];
        unsigned short nearEnd = 1 - FarEnd(tjs, pfp, start);
        float sep = PosSep2(pfp.XYZ[nearEnd], start);
        if(sep > close) continue;
        close = sep;
        pfpid = pid;
      }
      return pfpid;
    } // part
    return 0;
  } // PrimaryElectronPFPID

bool tca::MCParticleListUtils::PrimaryElectronStart	(	Point3_t &	start,
		Vector3_t &	dir,
		float &	energy
	)

Definition at line 1346 of file TCTruth.cxx.

References sim::ParticleList::EveId(), tca::TjStuff::MCPartList, part, cheat::ParticleInventoryService::ParticleList(), tca::SetMag(), and tca::TruthMatcher::tjs.

  {
    // returns the SCE corrected start position of a primary electron
    if(tjs.MCPartList.empty()) return false;
    art::ServiceHandle<cheat::ParticleInventoryService> pi_serv;
    for(unsigned int part = 0; part < tjs.MCPartList.size(); ++part) {
      auto& mcp = tjs.MCPartList[part];
      // require electron
      if(abs(mcp->PdgCode()) != 11) continue;
      int eveID = pi_serv->ParticleList().EveId(mcp->TrackId());
      if(mcp->TrackId() != eveID) continue;
      start = {{mcp->Vx(), mcp->Vy(), mcp->Vz()}};
      auto const* SCE = lar::providerFrom<spacecharge::SpaceChargeService>();
      geo::Vector_t posOffsets = SCE->GetPosOffsets({start[0], start[1], start[2]});
      posOffsets.SetX(-posOffsets.X());
      start[0] += posOffsets.X();
      start[1] += posOffsets.Y();
      start[2] += posOffsets.Z();
      dir = {{mcp->Px(), mcp->Py(), mcp->Pz()}};
      SetMag(dir, 1);
      energy = 1000 * (mcp->E() - mcp->Mass());
      return true;
    } // part
    return false;
  } // PrimaryElectronStart

int tca::MCParticleListUtils::PrimaryElectronTjID

(

CTP_t

inCTP

)

Definition at line 1437 of file TCTruth.cxx.

References sim::ParticleList::EveId(), tca::TjStuff::MCPartList, part, cheat::ParticleInventoryService::ParticleList(), and tca::TruthMatcher::tjs.

  {
    // returns the ID of a tj in inCTP that is closest to the start of a primary electron
    if(tjs.MCPartList.empty()) return 0;
    art::ServiceHandle<cheat::ParticleInventoryService> pi_serv;
    for(unsigned int part = 0; part < tjs.MCPartList.size(); ++part) {
      auto& mcp = tjs.MCPartList[part];
      // require electron
      if(abs(mcp->PdgCode()) != 11) continue;
      int eveID = pi_serv->ParticleList().EveId(mcp->TrackId());
      if(mcp->TrackId() != eveID) continue;
      return MCParticleStartTjID(part, inCTP);
    } // part
    return 0;
  } // PrimaryElectronTjID

Member Data Documentation

TjStuff& tca::MCParticleListUtils::tjs

Definition at line 86 of file TCTruth.h.

---

The documentation for this class was generated from the following files:

• larreco/v06_64_02/source/larreco/RecoAlg/TCAlg/TCTruth.h
• larreco/v06_64_02/source/larreco/RecoAlg/TCAlg/TCTruth.cxx