LArSoft  v06_85_00
Liquid Argon Software toolkit - http://larsoft.org/
TCTruth.cxx
Go to the documentation of this file.
3 
4 
7 
12 
16 
17 namespace tca {
18 
21  {
22  // Initialize the variables used to calculate Efficiency * Purity (aka EP) for matching to truth
23  EPCnts.fill(0);
24  TSums.fill(0.0);
25  EPTSums.fill(0.0);
26  TruVxCounts.fill(0);
27  nBadEP = 0;
28  } // Initialize
29 
32  {
33  // Matches reco hits to MC true tracks and puts the association into
34  // TCHit TruTrkID. This code is almost identical to the first part of MatchTruth.
35 
36  tjs.MCPartList.clear();
37 
38  if(tjs.MatchTruth[0] < 0) return;
39  if(tjs.fHits.empty()) return;
40 
43  // list of all true particles
44  sim::ParticleList const& plist = pi_serv->ParticleList();
45  if(plist.empty()) return;
46 
47  // MC Particles for the desired true particles
48  int sourcePtclTrackID = -1;
49 
50  // first find the source particle
51  simb::Origin_t sourceOrigin = simb::kUnknown;
52  for(sim::ParticleList::const_iterator ipart = plist.begin(); ipart != plist.end(); ++ipart) {
53  simb::MCParticle* mcp = (*ipart).second;
54  int trackID = mcp->TrackId();
55  art::Ptr<simb::MCTruth> theTruth = pi_serv->TrackIdToMCTruth_P(trackID);
56  if(tjs.MatchTruth[0] == 1) {
57  // Look for beam neutrino or single particle
58  if(theTruth->Origin() == simb::kBeamNeutrino) {
59  sourcePtclTrackID = trackID;
60  sourceOrigin = simb::kBeamNeutrino;
61  if(tjs.MatchTruth[1] > 0) {
62  Vector3_t dir {{mcp->Px(), mcp->Py(), mcp->Pz()}};
63  SetMag(dir, 1);
64  std::cout<<"Found beam neutrino sourcePtclTrackID "<<trackID<<" PDG code "<<mcp->PdgCode();
65  std::cout<<" Vx "<<std::fixed<<std::setprecision(1)<<mcp->Vx()<<" Vy "<<mcp->Vy()<<" Vz "<<mcp->Vz();
66  std::cout<<" dir "<<std::setprecision(3)<<dir[0]<<" "<<dir[1]<<" "<<dir[2]<<"\n";
67  }
68  break;
69  } // beam neutrino
70  if(theTruth->Origin() == simb::kSingleParticle) {
71  sourcePtclTrackID = trackID;
72  sourceOrigin = simb::kSingleParticle;
73  if(tjs.MatchTruth[1] > 0) {
74  Vector3_t dir {{mcp->Px(), mcp->Py(), mcp->Pz()}};
75  SetMag(dir, 1);
76  std::cout<<"Found single particle sourcePtclTrackID "<<trackID<<" PDG code "<<mcp->PdgCode();
77  std::cout<<" Vx "<<std::fixed<<std::setprecision(1)<<mcp->Vx()<<" Vy "<<mcp->Vy()<<" Vz "<<mcp->Vz();
78  std::cout<<" dir "<<std::setprecision(3)<<dir[0]<<" "<<dir[1]<<" "<<dir[2]<<"\n";
79  break;
80  }
81  } // single particle
82  } else if(tjs.MatchTruth[0] == 2 && theTruth->Origin() == simb::kCosmicRay) {
83  sourcePtclTrackID = trackID;
84  sourceOrigin = simb::kCosmicRay;
85  }
86  } // ipart
87 
88  if(sourcePtclTrackID == -1) {
89  if(tjs.MatchTruth[1] > 0) std::cout<<"MatchTrueHits: SourcePtcl not found\n";
90  return;
91  }
92 
93  if(tjs.MatchTruth[1] > 2) {
94  // print out a whole bunch of information
95  mf::LogVerbatim myprt("TC");
96  myprt<<"Displaying all neutrino origin MCParticles with T > 50 MeV\n";
97  myprt<<" trackID PDGCode Mother T(MeV) ________dir_______ Process\n";
98  for(sim::ParticleList::const_iterator ipart = plist.begin(); ipart != plist.end(); ++ipart) {
99  simb::MCParticle* mcp = (*ipart).second;
100  int trackID = mcp->TrackId();
101  art::Ptr<simb::MCTruth> theTruth = pi_serv->TrackIdToMCTruth_P(trackID);
102  if(theTruth->Origin() != simb::kBeamNeutrino) continue;
103  // Kinetic energy in MeV
104  int TMeV = 1000 * (mcp->E() - mcp->Mass());
105  if(TMeV < 50) continue;
106  myprt<<std::setw(8)<<trackID;
107  myprt<<std::setw(8)<<mcp->PdgCode();
108  myprt<<std::setw(8)<<mcp->Mother();
109  myprt<<std::setw(6)<<TMeV;
110  Point3_t pos;
111  pos[0] = mcp->Vx();
112  pos[1] = mcp->Vy();
113  pos[2] = mcp->Vz();
114  Vector3_t dir {{mcp->Px(), mcp->Py(), mcp->Pz()}};
115  SetMag(dir, 1);
116  myprt<<std::setprecision(2);
117  for(unsigned short xyz = 0; xyz < 3; ++xyz) myprt<<std::setw(6)<<dir[xyz];
118  myprt<<std::setw(22)<<mcp->Process();
119  myprt<<"\n";
120  } // ipart
121  } // big print
122 
123  // flag MCParticles to select for measuring performance.
124  std::vector<bool> select(plist.size(), false);
125  // ensure that we get all of the primary particles
126  unsigned int indx = 0;
127  for(sim::ParticleList::const_iterator ipart = plist.begin(); ipart != plist.end(); ++ipart) {
128  simb::MCParticle* mcp = (*ipart).second;
129  art::Ptr<simb::MCTruth> theTruth = pi_serv->TrackIdToMCTruth_P(mcp->TrackId());
130  if(theTruth->Origin() != sourceOrigin) continue;
131  select[indx] = true;
132  ++indx;
133  }
134 
135  // make a temp vector of hit -> geant trackID
136  std::vector<int> gtid(tjs.fHits.size(), 0);
137  // find hits that match to the source particle
138  for(unsigned int iht = 0; iht < tjs.fHits.size(); ++iht) {
139  std::vector<sim::TrackIDE> ides;
140  auto& tcHit = tjs.fHits[iht];
141  raw::ChannelID_t channel = tjs.geom->PlaneWireToChannel(tcHit.ArtPtr->WireID());
142  geo::PlaneID planeID = geo::PlaneID(tcHit.ArtPtr->WireID().Cryostat, tcHit.ArtPtr->WireID().TPC, tcHit.ArtPtr->WireID().Plane);
143  auto rhit = recob::Hit(channel,
144  tcHit.StartTick, tcHit.EndTick,
145  tcHit.PeakTime, tcHit.SigmaPeakTime,
146  tcHit.RMS,
147  tcHit.PeakAmplitude, tcHit.SigmaPeakAmp,
148  tcHit.Integral, tcHit.Integral, tcHit.SigmaIntegral,
149  tcHit.Multiplicity, tcHit.LocalIndex,
150  tcHit.GoodnessOfFit, tcHit.NDOF,
151  tjs.geom->View(channel),
152  tjs.geom->SignalType(planeID),
153  tcHit.ArtPtr->WireID());
154  try {
155  ides = bt_serv->HitToTrackIDEs(rhit);
156  }
157  catch(...) {}
158  if(ides.empty()) continue;
159  float energy = 0;
160  for(auto ide : ides) energy += ide.energy;
161  if(energy == 0) continue;
162  // require 1/2 of the energy be due to one MC particle
163  energy /= 2;
164  int hitTruTrkID = 0;
165  for(auto ide : ides) {
166  if(ide.energy > energy) {
167  hitTruTrkID = ide.trackID;
168  break;
169  }
170  } // ide
171  if(hitTruTrkID == 0) continue;
172  // ensure it has the correct source
173  art::Ptr<simb::MCTruth> theTruth = pi_serv->TrackIdToMCTruth_P(hitTruTrkID);
174  if(theTruth->Origin() != sourceOrigin) continue;
175  unsigned int indx = 0;
176  for(sim::ParticleList::const_iterator ipart = plist.begin(); ipart != plist.end(); ++ipart) {
177  ++indx;
178  simb::MCParticle* mcp = (*ipart).second;
179  if(mcp->TrackId() != hitTruTrkID) continue;
180  select[indx - 1] = true;
181  gtid[iht] = mcp->TrackId();
182  // find the eve particle and select it as well
183  const simb::MCParticle* momMCP = pi_serv->TrackIdToMotherParticle_P(mcp->TrackId());
184  if(!momMCP) continue;
185  if(momMCP->TrackId() == mcp->TrackId()) break;
186  unsigned int mindx = 0;
187  for(sim::ParticleList::const_iterator mpart = plist.begin(); mpart != plist.end(); ++mpart) {
188  simb::MCParticle* mmcp = (*mpart).second;
189  if(mmcp->TrackId() == momMCP->TrackId()) {
190  select[mindx] = true;
191  break;
192  }
193  ++mindx;
194  } // mpart
195  break;
196  } // ipart
197  } // iht
198 
199  // save the selected MCParticles in tjs
200  indx = 0;
201  for(sim::ParticleList::const_iterator ipart = plist.begin(); ipart != plist.end(); ++ipart) {
202  if(select[indx]) {
203  simb::MCParticle* mcp = (*ipart).second;
204  tjs.MCPartList.push_back(mcp);
205  }
206  ++indx;
207  } // ipart
208 
209  if(tjs.MCPartList.size() > UINT_MAX) {
210  std::cout<<"MatchTrueHits: Crazy large number of MCParticles "<<tjs.MCPartList.size()<<". Ignoring this event\n";
211  tjs.MCPartList.clear();
212  return;
213  }
214 
215  // define MCPartListIndex for the hits
216  unsigned int nMatch = 0;
217  for(unsigned int iht = 0; iht < tjs.fHits.size(); ++iht) {
218  if(gtid[iht] == 0) continue;
219  auto& hit = tjs.fHits[iht];
220  for(unsigned int indx = 0; indx < tjs.MCPartList.size(); ++indx) {
221  auto& mcp = tjs.MCPartList[indx];
222  if(mcp->TrackId() != gtid[iht]) continue;
223  hit.MCPartListIndex = indx;
224  ++nMatch;
225  break;
226  } // indx
227  } // iht
228 
229 
230  std::cout<<"MatchTrueHits: MCPartList size "<<tjs.MCPartList.size()<<" nMatched hits "<<nMatch<<"\n";
231 
232  } // MatchTrueHits
233 
236  {
237  // study characteristics of shower parent pfps. This code is adapted from TCShower FindParent
238  if(tjs.pfps.empty()) return;
239  if(tjs.MCPartList.empty()) return;
240 
241  // Look for truth pfp primary electron
242  Point3_t primVx {{-666.0, -666.0, -666.0}};
243  // the primary should be the first one in the list as selected in GetHitCollection
244  auto& primMCP = tjs.MCPartList[0];
245  primVx[0] = primMCP->Vx();
246  primVx[1] = primMCP->Vy();
247  primVx[2] = primMCP->Vz();
248  geo::Vector_t posOffsets;
249  auto const* SCE = lar::providerFrom<spacecharge::SpaceChargeService>();
250  posOffsets = SCE->GetPosOffsets({primVx[0], primVx[1], primVx[2]});
251  posOffsets.SetX(-posOffsets.X());
252  primVx[0] += posOffsets.X();
253  primVx[1] += posOffsets.Y();
254  primVx[2] += posOffsets.Z();
255  geo::TPCID inTPCID = tjs.TPCID;
256  if(!InsideTPC(tjs, primVx, inTPCID)) return;
257 
258  std::string fcnLabel = "SSP";
259  // Create a truth shower for each primary electron
261  MCParticleListUtils mcpu{tjs};
262  for(unsigned int part = 0; part < tjs.MCPartList.size(); ++part) {
263  auto& mcp = tjs.MCPartList[part];
264  // require electron or photon
265  if(abs(mcp->PdgCode()) != 11 && abs(mcp->PdgCode()) != 111) continue;
266  int eveID = pi_serv->ParticleList().EveId(mcp->TrackId());
267  // require that it is primary
268  if(mcp->TrackId() != eveID) continue;
269  int truPFP = 0;
270  auto ss3 = mcpu.MakeCheatShower(part, primVx, truPFP);
271  if(ss3.ID == 0) continue;
272  if(truPFP == 0) continue;
273  if(!StoreShower(fcnLabel, tjs, ss3)) {
274  std::cout<<"Failed to store 3S"<<ss3.ID<<"\n";
275  break;
276  } // store failed
277  // now fill the TTree
278  float ss3Energy = ShowerEnergy(ss3);
279  for(auto& pfp : tjs.pfps) {
280  if(pfp.TPCID != ss3.TPCID) continue;
281  // ignore neutrinos
282  if(pfp.PDGCode == 12 || pfp.PDGCode == 14) continue;
283  // ignore shower pfps
284  if(pfp.PDGCode == 1111) continue;
285  float pfpEnergy = 0;
286  float minEnergy = 1E6;
287  for(auto tid : pfp.TjIDs) {
288  auto& tj = tjs.allTraj[tid - 1];
289  float energy = ChgToMeV(tj.TotChg);
290  pfpEnergy += energy;
291  if(energy < minEnergy) minEnergy = energy;
292  }
293  pfpEnergy -= minEnergy;
294  pfpEnergy /= (float)(pfp.TjIDs.size() - 1);
295  // find the end that is farthest away
296  unsigned short pEnd = FarEnd(tjs, pfp, ss3.ChgPos);
297  auto pToS = PointDirection(pfp.XYZ[pEnd], ss3.ChgPos);
298  // take the absolute value in case the shower direction isn't well known
299  float costh1 = std::abs(DotProd(pToS, ss3.Dir));
300  float costh2 = DotProd(pToS, pfp.Dir[pEnd]);
301  // distance^2 between the pfp end and the shower start, charge center, and shower end
302  float distToStart2 = PosSep2(pfp.XYZ[pEnd], ss3.Start);
303  float distToChgPos2 = PosSep2(pfp.XYZ[pEnd], ss3.ChgPos);
304  float distToEnd2 = PosSep2(pfp.XYZ[pEnd], ss3.End);
305 // mf::LogVerbatim("TC")<<" 3S"<<ss3.ID<<" P"<<pfp.ID<<"_"<<pEnd<<" distToStart "<<sqrt(distToStart2)<<" distToChgPos "<<sqrt(distToChgPos2)<<" distToEnd "<<sqrt(distToEnd2);
306  // find the end of the shower closest to the pfp
307  unsigned short shEnd = 0;
308  if(distToEnd2 < distToStart2) shEnd = 1;
309  if(shEnd == 0 && distToChgPos2 < distToStart2) continue;
310  if(shEnd == 1 && distToChgPos2 < distToEnd2) continue;
311 // mf::LogVerbatim("TC")<<" 3S"<<ss3.ID<<"_"<<shEnd<<" P"<<pfp.ID<<"_"<<pEnd<<" costh1 "<<costh1;
312  Point2_t alongTrans;
313  // find the longitudinal and transverse components of the pfp start point relative to the
314  // shower center
315  FindAlongTrans(ss3.ChgPos, ss3.Dir, pfp.XYZ[pEnd], alongTrans);
316 // mf::LogVerbatim("TC")<<" alongTrans "<<alongTrans[0]<<" "<<alongTrans[1];
317  hist.fSep = sqrt(distToChgPos2);
318  hist.fShEnergy = ss3Energy;
319  hist.fPfpEnergy = pfpEnergy;
320  hist.fPfpLen = PosSep(pfp.XYZ[0], pfp.XYZ[1]);
321  hist.fMCSMom = MCSMom(tjs, pfp.TjIDs);
322  hist.fDang1 = acos(costh1);
323  hist.fDang2 = acos(costh2);
324  hist.fChgFrac = 0;
325  float chgFrac = 0;
326  float totSep = 0;
327  // find the charge fraction btw the pfp start and the point that is
328  // half the distance to the charge center in each plane
329  for(unsigned short plane = 0; plane < tjs.NumPlanes; ++plane) {
330  CTP_t inCTP = EncodeCTP(ss3.TPCID.Cryostat, ss3.TPCID.TPC, plane);
331  int ssid = 0;
332  for(auto cid : ss3.CotIDs) {
333  auto& ss = tjs.cots[cid - 1];
334  if(ss.CTP != inCTP) continue;
335  ssid = ss.ID;
336  break;
337  } // cid
338  if(ssid == 0) continue;
339  auto tpFrom = MakeBareTP(tjs, pfp.XYZ[pEnd], pToS, inCTP);
340  auto& ss = tjs.cots[ssid - 1];
341  auto& stp1 = tjs.allTraj[ss.ShowerTjID - 1].Pts[1];
342  float sep = PosSep(tpFrom.Pos, stp1.Pos);
343  float toPos = tpFrom.Pos[0] + 0.5 * tpFrom.Dir[0] * sep;
344  float cf = ChgFracBetween(tjs, tpFrom, toPos, false);
345  // weight by the separation in the plane
346  totSep += sep;
347  chgFrac += sep * cf;
348  } // plane
349  if(totSep > 0) hist.fChgFrac = chgFrac / totSep;
350  hist.fAlong = alongTrans[0];
351  hist.fTrans = alongTrans[1];
352  hist.fInShwrProb = InShowerProbLong(ss3Energy, -hist.fSep);
353  bool isBad = (hist.fDang1 > 2 || hist.fChgFrac < 0.5 || hist.fInShwrProb < 0.05);
354  if(pfp.ID == truPFP && isBad) {
355  mf::LogVerbatim myprt("TC");
356  myprt<<"SSP: 3S"<<ss3.ID<<" shEnergy "<<(int)ss3Energy<<" P"<<pfp.ID<<" pfpEnergy "<<(int)pfpEnergy;
357  myprt<<" MCSMom "<<hist.fMCSMom<<" len "<<hist.fPfpLen;
358  myprt<<" Dang1 "<<hist.fDang1<<" Dang2 "<<hist.fDang2<<" chgFrac "<<hist.fChgFrac;
359  myprt<<" fInShwrProb "<<hist.fInShwrProb;
360  myprt<<" EventsProcessed "<<tjs.EventsProcessed;
361  }
362  if(pfp.ID == truPFP) {
363  hist.fShowerParentSig->Fill();
364  } else {
365  hist.fShowerParentBkg->Fill();
366  }
367  } // pfp
368  } // part
369 // PrintShowers(fcnLabel, tjs);
370 // Print2DShowers(fcnLabel, tjs, USHRT_MAX, false);
371  // kill the cheat showers
372  for(auto& ss3 : tjs.showers) {
373  if(ss3.ID == 0) continue;
374  if(!ss3.Cheat) continue;
375  for(auto cid : ss3.CotIDs) {
376  auto& ss = tjs.cots[cid - 1];
377  ss.ID = 0;
378  auto& stj = tjs.allTraj[ss.ShowerTjID - 1];
379  stj.AlgMod[kKilled] = true;
380  } // cid
381  ss3.ID = 0;
382  } // ss3
383  } // StudyShowerParents
384 
387  {
388  // study tjs matched to electrons to develop an electron tag
389 
390 // float likely;
391 // bool flipDirection;
392  for(auto& tj : tjs.allTraj) {
393  if(tj.AlgMod[kKilled]) continue;
394  if(tj.MCPartListIndex == UINT_MAX) continue;
395  unsigned short npts = tj.EndPt[1] - tj.EndPt[0] + 1;
396  if(npts < 10) continue;
397 // PrimaryElectronLikelihood(tjs, tj, likely, flipDirection, true);
398  auto& mcp = tjs.MCPartList[tj.MCPartListIndex];
399  unsigned short pdgIndex = PDGCodeIndex(tjs, mcp->PdgCode());
400  if(pdgIndex > 4) continue;
401  unsigned short midpt = 0.5 * (tj.EndPt[0] + tj.EndPt[1]);
402  float mom1 = MCSMom(tjs, tj, tj.EndPt[0], midpt);
403  float mom2 = MCSMom(tjs, tj, midpt, tj.EndPt[1]);
404  float asym = std::abs(mom1 - mom2) / (mom1 + mom2);
405  hist.fChgRMS[pdgIndex]->Fill(tj.ChgRMS);
406  hist.fMomAsym[pdgIndex]->Fill(asym);
407  float elike = asym * tj.ChgRMS;
408  hist.fElectronLike[pdgIndex]->Fill(elike);
409  float len = PosSep(tj.Pts[tj.EndPt[0]].Pos, tj.Pts[tj.EndPt[1]].Pos);
410  hist.fElectronLike_Len[pdgIndex]->Fill(len, elike);
411  } // tj
412  } // StudyElectrons
413 
416  {
419  sim::ParticleList const& plist = pi_serv->ParticleList();
420 
421  unsigned short cnt = 0;
422  for(sim::ParticleList::const_iterator ipart = plist.begin(); ipart != plist.end(); ++ipart) {
423  ++cnt;
424  const simb::MCParticle* p = (*ipart).second;
425  int pdg = abs(p->PdgCode());
426  if(cnt == 1 && pdg != 111) {
427  std::cout<<"First MC particle isn't a pizero\n";
428  return;
429  }
430  if(cnt == 1) continue;
431  if(pdg != 22) break;
432  double photE = 1000 * p->E();
433  if(photE < 50) continue;
434 // Point3_t photStart {p->Vx(), p->Vy(), p->Vz()};
435  Vector3_t photDir {{p->Px(), p->Py(), p->Pz()}};
436  SetMag(photDir, 1);
437  // sum up the charge for all hits that are daughters of this photon
438  std::vector<float> chgSum(3);
439  for(unsigned int iht = 0; iht < tjs.fHits.size(); ++iht) {
440  auto mcpIndex = tjs.fHits[iht].MCPartListIndex;
441  if(mcpIndex == UINT_MAX) continue;
442  auto& mcp = tjs.MCPartList[mcpIndex];
443  int eveID = pi_serv->ParticleList().EveId(mcp->TrackId());
444  if(eveID != p->TrackId()) continue;
445  unsigned short plane = tjs.fHits[iht].ArtPtr->WireID().Plane;
446  chgSum[plane] += tjs.fHits[iht].Integral;
447  } // iht
448  for(unsigned short plane = 0; plane < 3; ++plane) {
449  if(chgSum[plane] == 0) continue;
450  float chgCal = photE / chgSum[plane];
451  hist.fChgToMeV[plane]->Fill(chgCal);
452  hist.fChgToMeV_Etru->Fill(photE, chgCal);
453  } // plane
454  } // ipart
455 
456  } // StudyPiZeros
457 
459  void TruthMatcher::MatchTruth(const HistStuff& hist, bool fStudyMode)
460  {
461  // The hits have already been matched to the truth in MatchTrueHits. Here we match reconstructed objects
462  // to the truth-matched hits to measure performance
463 
464  if(tjs.MatchTruth[0] < 0) return;
465  if(tjs.MCPartList.empty()) return;
466  for(auto& pfp : tjs.pfps) pfp.EffPur = 0;
467 
469  geo::Vector_t posOffsets;
470  auto const* SCE = lar::providerFrom<spacecharge::SpaceChargeService>();
471 // if(!SCE->EnableSimSpatialSCE()) std::cout<<">>>> EnableSimSpatialSCE isn't active...\n";
472 
473  // these are only applicable to neutrinos
474  bool neutrinoVxReconstructable = false;
475  bool vxReconstructedNearNuVx = false;
476  bool neutrinoPFPCorrect = false;
477  // Locate the primary vertex
478  Point3_t primVx {{-666.0, -666.0, -666.0}};
479  // the primary should be the first one in the list as selected in GetHitCollection
480  auto& primMCP = tjs.MCPartList[0];
481  primVx[0] = primMCP->Vx();
482  primVx[1] = primMCP->Vy();
483  primVx[2] = primMCP->Vz();
484  posOffsets = SCE->GetPosOffsets({primVx[0], primVx[1], primVx[2]});
485  posOffsets.SetX(-posOffsets.X());
486  primVx[0] += posOffsets.X();
487  primVx[1] += posOffsets.Y();
488  primVx[2] += posOffsets.Z();
489  if(tjs.MatchTruth[1] > 1) std::cout<<"Prim PDG code "<<primMCP->PdgCode()<<" primVx "<<std::fixed<<std::setprecision(1)<<primVx[0]<<" "<<primVx[1]<<" "<<primVx[2]<<"\n";
490  geo::TPCID inTPCID = tjs.TPCID;
491  if(!InsideTPC(tjs, primVx, inTPCID)) {
492  if(tjs.MatchTruth[1] > 0) std::cout<<"Found a primary particle but it is not inside any TPC\n";
493  return;
494  }
495  neutrinoVxReconstructable = true;
496 
497  // Look for the MC truth process that should be considered (beam neutrino,
498  // single particle, cosmic rays), then form a list of selected MCParticles
499  // that will be used to measure performance. Feb 16: Changed GetHitCollection to only
500  // store the MCParticles for the desired MCTruth collection
501  std::vector<unsigned int> mcpSelect;
502  // vector of reconstructable primary particles
503  std::vector<unsigned int> primMCPs;
504  for(unsigned int part = 0; part < tjs.MCPartList.size(); ++part) {
505  auto& mcp = tjs.MCPartList[part];
506  // require it is charged
507  int pdg = abs(mcp->PdgCode());
508  bool isCharged = (pdg == 11) || (pdg == 13) || (pdg == 211) || (pdg == 321) || (pdg == 2212);
509  if(!isCharged) continue;
510  // require that it can be reconstructed in 3D
511  if(!CanReconstruct(part, 3, inTPCID)) continue;
512  mcpSelect.push_back(part);
513  // Now require MCParticle primaries
514  if(mcp->Mother() != 0) continue;
515  // add to the list of primaries
516  primMCPs.push_back(part);
517  } // part
518  if(mcpSelect.empty()) return;
519  tjs.SelectEvent = true;
520 // mf::LogVerbatim("TC")<<"SelectEvent "<<tjs.Run<<" "<<tjs.SubRun<<" "<<tjs.Event;
521 
522  if(neutrinoVxReconstructable) ++TruVxCounts[0];
523 
524  // Form a list of mother-daughter pairs that should be considered as a single particle
525  std::vector<std::pair<unsigned int, unsigned int>> moda;
526  for(unsigned int part = 0; part < mcpSelect.size(); ++part) {
527  auto& mcp = tjs.MCPartList[mcpSelect[part]];
528  if(mcp->NumberDaughters() == 0) continue;
529  unsigned int ndtr = 0;
530  unsigned int dtrSelIndex = 0;
531  for(int idtr = 0; idtr < mcp->NumberDaughters(); ++idtr) {
532  int dtrTrackId = mcp->Daughter(idtr);
533  // ignore it if it's not in the list
534  bool ignore = true;
535  for(unsigned short dpart = part + 1; dpart < mcpSelect.size(); ++dpart) {
536  auto& dmcp = tjs.MCPartList[mcpSelect[dpart]];
537  if(dmcp->TrackId() == dtrTrackId) {
538  dtrSelIndex = dpart;
539  ignore = false;
540  }
541  } // dpart
542  if(ignore) continue;
543  ++ndtr;
544  } // idtr
545  // require only one daughter
546  if(ndtr != 1) continue;
547  // require that the daughter have the same PDG code
548  auto& dtr = tjs.MCPartList[mcpSelect[dtrSelIndex]];
549  if(dtr->PdgCode() != mcp->PdgCode()) continue;
550  if(tjs.MatchTruth[1] > 1) mf::LogVerbatim("TC")<<"Daughter MCP "<<dtrSelIndex<<" -> Mother MCP "<<part;
551  moda.push_back(std::make_pair(mcpSelect[part], mcpSelect[dtrSelIndex]));
552  // delete the daughter entry
553  mcpSelect.erase(mcpSelect.begin() + dtrSelIndex);
554  } // part
555 
556  // NOTE: The PutMCPHitsInVector function will return an incomplete set of hits
557  // matched to a MCParticle if it is called before mother-daughter assocations are corrected below.
558  if(!moda.empty()) {
559  // over-write the daughter hit -> MCParticle association with the mother.
560  // Note that mother-daughter pairs are ordered by increasing generation. Reverse
561  // moda so that the grand-daughters come first. Grand-daughters will then be
562  // over-written by daughters in the previous generation
563  if(moda.size() > 1) std::reverse(moda.begin(), moda.end());
564  for(auto& hit : tjs.fHits) {
565  // see if this is a daughter
566  for(auto& md : moda) if(md.second == hit.MCPartListIndex) hit.MCPartListIndex = md.first;
567  } // hit
568  } // !moda.empty
569 
570  // True histograms
571  for(unsigned int part = 0; part < mcpSelect.size(); ++part) {
572  auto& mcp = tjs.MCPartList[mcpSelect[part]];
573  unsigned short pdgIndex = PDGCodeIndex(tjs, mcp->PdgCode());
574  if(pdgIndex > 4) continue;
575  float TMeV = 1000 * (mcp->E() - mcp->Mass());
576  hist.fTruT[pdgIndex]->Fill(TMeV);
577  } // part
578 
579  // Match tjs to MC particles. Declare it a match to the MC particle that has the most
580  // hits in the tj
581  std::vector<std::vector<int>> tjid(tjs.NumPlanes);
582  std::vector<std::vector<unsigned short>> ntht(tjs.NumPlanes);
583  for(unsigned short plane = 0; plane < tjs.NumPlanes; ++plane) {
584  tjid[plane].resize(tjs.MCPartList.size());
585  ntht[plane].resize(tjs.MCPartList.size());
586  }
587 
588  for(auto& tj : tjs.allTraj) {
589  if(tj.AlgMod[kKilled]) continue;
590  geo::PlaneID planeID = DecodeCTP(tj.CTP);
591  if(planeID.Cryostat != inTPCID.Cryostat) continue;
592  if(planeID.TPC != inTPCID.TPC) continue;
593  tj.MCPartListIndex = UINT_MAX;
594  std::vector<unsigned int> mcpIndex, cnt;
595  auto tjhits = PutTrajHitsInVector(tj, kUsedHits);
596  for(auto iht : tjhits) {
597  auto& hit = tjs.fHits[iht];
598  if(hit.MCPartListIndex > tjs.MCPartList.size() - 1) continue;
599  // ignore Tj hits that aren't in mcpSelect
600  if(std::find(mcpSelect.begin(), mcpSelect.end(), hit.MCPartListIndex) == mcpSelect.end()) continue;
601  unsigned int indx = 0;
602  for(indx = 0; indx < mcpIndex.size(); ++indx) if(hit.MCPartListIndex == mcpIndex[indx]) break;
603  if(indx == mcpIndex.size()) {
604  mcpIndex.push_back(hit.MCPartListIndex);
605  cnt.push_back(1);
606  } else {
607  ++cnt[indx];
608  }
609  } // iht
610  unsigned short maxcnt = 0;
611  unsigned int tmpIndex = UINT_MAX;
612  for(unsigned short ii = 0; ii < cnt.size(); ++ii) {
613  if(cnt[ii] > maxcnt) {
614  maxcnt = cnt[ii];
615  tmpIndex = mcpIndex[ii];
616  }
617  } // ii
618  if(tmpIndex == UINT_MAX) continue;
619  // calculate efficiency and purity.
620  // Put the truth-matched hits into a vector for counting
621  auto mcpHits = PutMCPHitsInVector(tmpIndex, tj.CTP);
622  if(mcpHits.size() < 3) continue;
623  float eff = (float)maxcnt / (float)mcpHits.size();
624  float pur = (float)maxcnt / (float)tjhits.size();
625  float ep = eff * pur;
626  // see if there is a previous match with more truth-matched hits
627  unsigned short plane = DecodeCTP(tj.CTP).Plane;
628  if(tjid[plane][tmpIndex] > 0) {
629  if(maxcnt > ntht[plane][tmpIndex]) {
630  // clear the old one
631  auto& mtj = tjs.allTraj[tjid[plane][tmpIndex] - 1];
632  mtj.EffPur = 0;
633  mtj.MCPartListIndex = UINT_MAX;
634  // use the new one
635  tj.EffPur = ep;
636  tj.MCPartListIndex = tmpIndex;
637  tjid[plane][tmpIndex] = tj.ID;
638  ntht[plane][tmpIndex] = maxcnt;
639  } else {
640  // use the old one
641  continue;
642  }
643  } else {
644  // no previous match
645  tj.EffPur = eff * pur;
646  tj.MCPartListIndex = tmpIndex;
647  tjid[plane][tmpIndex] = tj.ID;
648  ntht[plane][tmpIndex] = maxcnt;
649  }
650  } // tj
651 
652  if(neutrinoVxReconstructable) {
653  // find the vertex closest to the true primary vertex
654  float best = 1;
655  unsigned short vx3ID = 0;
656  if(!tjs.pfps.empty()) {
657  auto& pfp = tjs.pfps[0];
658  if((pfp.PDGCode == 14 || pfp.PDGCode == 12) && pfp.Vx3ID[0] > 0) {
659  // Found a neutrino pfp with a start vertex
660  auto& vx3 = tjs.vtx3[pfp.Vx3ID[0] - 1];
661  // check the proximity to the true vertex
662  // BUG the double brace syntax is required to work around clang bug 21629
663  // (https://bugs.llvm.org/show_bug.cgi?id=21629)
664  Point3_t vpos = {{vx3.X, vx3.Y, vx3.Z}};
665  if(PosSep(vpos, primVx) < 1) neutrinoPFPCorrect = true;
666  } // neutrino pfp
667  } // PFParticles exist
668  for(auto& vx3 : tjs.vtx3) {
669  if(vx3.ID == 0) continue;
670  if(vx3.TPCID != inTPCID) continue;
671  // BUG the double brace syntax is required to work around clang bug 21629
672  // (https://bugs.llvm.org/show_bug.cgi?id=21629)
673  Point3_t vpos = {{vx3.X, vx3.Y, vx3.Z}};
674  float sep = PosSep(vpos, primVx);
675  if(sep < best) {
676  best = sep;
677  vx3ID = vx3.ID;
678  }
679  } // vx3
680  if(vx3ID > 0) vxReconstructedNearNuVx = true;
681  } // neutrinoVxReconstructable
682 
683  if(tjs.MatchTruth[1] > 0) {
684  // print out
685  mf::LogVerbatim myprt("TC");
686  myprt<<"Number of primary particles "<<primMCPs.size()<<" Vtx";
687  for(unsigned short xyz = 0; xyz < 3; ++xyz) myprt<<" "<<std::fixed<<std::setprecision(1)<<primVx[xyz];
688  myprt<<" nuVx Reconstructable? "<<neutrinoVxReconstructable<<" vx near nuVx? "<<vxReconstructedNearNuVx;
689  myprt<<" neutrinoPFPCorrect? "<<neutrinoPFPCorrect<<"\n";
690  myprt<<"MCPIndx TrackId PDG eveID KE Len _______Dir_______ ____ProjInPln___ Process StartHit-EndHit_nTruHits";
691  for(unsigned int ipart = 0; ipart < tjs.MCPartList.size(); ++ipart) {
692  bool doPrt = (std::find(mcpSelect.begin(), mcpSelect.end(), ipart) != mcpSelect.end());
693  auto& mcp = tjs.MCPartList[ipart];
694  // also print if this is a pizero or decay photon > 30 MeV
695  if(mcp->PdgCode() == 111) doPrt = true;
696  // Kinetic energy in MeV
697  int TMeV = 1000 * (mcp->E() - mcp->Mass());
698  if(mcp->PdgCode() == 22 && mcp->Process() == "Decay" && TMeV > 30) doPrt = true;
699  if(!doPrt) continue;
700  myprt<<"\n";
701  myprt<<std::setw(7)<<ipart;
702  myprt<<std::setw(8)<<mcp->TrackId();
703  myprt<<std::setw(6)<<mcp->PdgCode();
704  myprt<<std::setw(8)<<pi_serv->ParticleList().EveId(mcp->TrackId());
705  myprt<<std::setw(6)<<TMeV;
706  Point3_t start {{mcp->Vx(), mcp->Vy(), mcp->Vz()}};
707  posOffsets = SCE->GetPosOffsets({start[0], start[1], start[2]});
708  posOffsets.SetX(-posOffsets.X());
709  start[0] += posOffsets.X();
710  start[1] += posOffsets.Y();
711  start[2] += posOffsets.Z();
712  Point3_t end {{mcp->EndX(), mcp->EndY(), mcp->EndZ()}};
713  posOffsets = SCE->GetPosOffsets({end[0], end[1], end[2]});
714  posOffsets.SetX(-posOffsets.X());
715  end[0] += posOffsets.X();
716  end[1] += posOffsets.Y();
717  end[2] += posOffsets.Z();
718  myprt<<std::setw(7)<<std::setprecision(1)<<PosSep(start, end);
719  Vector3_t dir {{mcp->Px(), mcp->Py(), mcp->Pz()}};
720  SetMag(dir, 1);
721  for(unsigned short xyz = 0; xyz < 3; ++xyz) myprt<<std::setw(6)<<std::setprecision(2)<<dir[xyz];
722  std::vector<float> startWire(tjs.NumPlanes);
723  for(unsigned short plane = 0; plane < tjs.NumPlanes; ++plane) {
724  CTP_t inCTP = plane;
725  auto tp = MakeBareTP(tjs, start, dir, inCTP);
726  myprt<<std::setw(6)<<tp.Delta;
727  startWire[plane] = tp.Pos[0];
728  } // plane
729  myprt<<std::setw(22)<<mcp->Process();
730  // print the extent of the particle in each TPC plane
731  for(const geo::TPCID& tpcid : tjs.geom->IterateTPCIDs()) {
732  geo::TPCGeo const& TPC = tjs.geom->TPC(tpcid);
733  for(unsigned short plane = 0; plane < TPC.Nplanes(); ++plane) {
734  CTP_t inCTP = EncodeCTP(tpcid.Cryostat, tpcid.TPC, plane);
735  auto mcpHits = PutMCPHitsInVector(ipart, inCTP);
736  if(mcpHits.empty()) continue;
737  // get the direction correct
738  float fwire = tjs.fHits[mcpHits[0]].ArtPtr->WireID().Wire;
739  float lwire = tjs.fHits[mcpHits[mcpHits.size() - 1]].ArtPtr->WireID().Wire;
740  if(std::abs(fwire - startWire[plane]) < std::abs(lwire - startWire[plane])) {
741  myprt<<" "<<PrintHitShort(tjs.fHits[mcpHits[0]])<<"-"<<PrintHitShort(tjs.fHits[mcpHits[mcpHits.size() - 1]]);
742  } else {
743  myprt<<" "<<PrintHitShort(tjs.fHits[mcpHits[mcpHits.size() - 1]])<<"-"<<PrintHitShort(tjs.fHits[mcpHits[0]]);
744  }
745  myprt<<"_"<<mcpHits.size();
746  } // plane
747  } // tpcid
748  } // ipart
749  } // tjs.MatchTruth[1] > 0
750 
751  // Match Tjs and PFParticles and accumulate statistics
752  MatchAndSum(hist, mcpSelect, inTPCID);
753 
754 /* This needs work...
755  // match 2D vertices (crudely)
756  for(auto& tj : tjs.allTraj) {
757  // obsolete vertex
758  if(tj.AlgMod[kKilled]) continue;
759  // require a truth match
760  if(tj.MCPartListIndex == UINT_MAX) continue;
761  if (tj.MCPartListIndex < tjs.MCPartList.size()) {
762  // ignore electrons unless it is a primary electron
763  auto& mcp = tjs.MCPartList[tj.MCPartListIndex];
764  int pdg = abs(mcp->PdgCode());
765  if(pdg == 11 && mcp->Mother() != 0) continue;
766  }
767  for(unsigned short end = 0; end < 2; ++end) {
768  if(tj.VtxID[end] == 0) continue;
769  VtxStore& vx2 = tjs.vtx[tj.VtxID[end]-1];
770  vx2.Stat[kVtxTruMatch] = true;
771  } // end
772  } // tj
773 */
774 
775  // Tj histograms
776  for(auto& tj : tjs.allTraj) {
777  if(tj.AlgMod[kKilled]) continue;
778  if(tj.MCPartListIndex == UINT_MAX) continue;
779  auto& mcp = tjs.MCPartList[tj.MCPartListIndex];
780  short truIndex = PDGCodeIndex(tjs, mcp->PdgCode());
781  if(truIndex == SHRT_MAX) continue;
782  float frac = 0;
783  float cnt = 0;
784  for(unsigned short ipt = tj.EndPt[0]; ipt <= tj.EndPt[1]; ++ipt) {
785  auto& tp = tj.Pts[ipt];
786  if(tp.Environment[kEnvNearTj]) ++frac;
787  ++cnt;
788  } // ipt
789  if(cnt > 0) frac /= cnt;
790  float TMeV = 1000 * (mcp->E() - mcp->Mass());
791  hist.fNearTj[truIndex]->Fill(TMeV, frac);
792  } // tj
793 
794  // match a 3D vertex to the primary vertex
795 /*
796  unsigned short vx3RecoPrmary = 0;
797  float close = 1E6;
798  for(auto& vx3 : tjs.vtx3) {
799  if(vx3.ID == 0) continue;
800  // ignore pfp start vertices
801  if(vx3.Wire == -2) continue;
802  Point3_t vxPos {{vx3.X, vx3.Y, vx3.Z}};
803  float sep = PosSep(vxPos, primVx);
804  if(sep < close) {
805  close = sep;
806  vx3RecoPrmary = vx3.ID;
807  }
808  } // vx3
809  if(vx3RecoPrmary > 0) {
810  auto& vx3 = tjs.vtx3[vx3RecoPrmary - 1];
811  std::cout<<"vx3RecoPrmary 3"<<vx3.ID<<" close "<<std::setprecision(1)<<close<<" deltas";
812  Point3_t vxPos {{vx3.X, vx3.Y, vx3.Z}};
813  for(unsigned short xyz = 0; xyz < 3; ++xyz) std::cout<<" "<<vxPos[xyz] - primVx[xyz];
814  std::cout<<"\n";
815  }
816 */
817 
818  // PFParticle histograms
819  constexpr double twopi = 2 * M_PI;
820  std::array<int, 5> recoCodeList = {{0, 11, 13, 211, 2212}};
821  for(auto& pfp : tjs.pfps) {
822  if(pfp.ID == 0) continue;
823  // ignore showers
824  if(pfp.PDGCode == 1111) continue;
825  // require match to MC
826  if(pfp.MCPartListIndex == UINT_MAX) continue;
827  auto& mcp = tjs.MCPartList[pfp.MCPartListIndex];
828  short truIndex = PDGCodeIndex(tjs, mcp->PdgCode());
829  if(truIndex == SHRT_MAX) continue;
830  short recIndex = 0;
831  for(recIndex = 0; recIndex < 5; ++recIndex) if(pfp.PDGCode == recoCodeList[recIndex]) break;
832  if(recIndex > 4) continue;
833  hist.PDGCode_reco_true->Fill((float)truIndex, (float)recIndex);
834  // get the true start position and shift it by the SCE offset
835  Point3_t truStart {{mcp->Vx(), mcp->Vy(), mcp->Vz()}};
836  Point3_t truEnd {{mcp->EndX(), mcp->EndY(), mcp->EndZ()}};
837  float truLen = PosSep(truStart, truEnd);
838  if(truLen < 2) continue;
839  posOffsets = SCE->GetPosOffsets({truStart[0], truStart[1], truStart[2]});
840  posOffsets.SetX(-posOffsets.X());
841  truStart[0] += posOffsets.X();
842  truStart[1] += posOffsets.Y();
843  truStart[2] += posOffsets.Z();
844  auto recDir = pfp.Dir[0];
845  short startEnd = 0;
846  if(PosSep(pfp.XYZ[1], truStart) < PosSep(pfp.XYZ[0], truStart)) {
847  startEnd = 1;
848  for(unsigned short xyz = 0; xyz < 3; ++xyz) recDir[xyz] *= -1;
849  }
850  hist.fPFPStartEnd->Fill((float)startEnd);
851 /*
852  if(std::abs(pfp.XYZ[startEnd][1] - truStart[1]) < 2 && std::abs(pfp.XYZ[startEnd][2] - truStart[2]) < 2) {
853  std::cout<<"MatchTruth P"<<pfp.ID<<" MCP "<<mcp->TrackId();
854  for(unsigned short xyz = 0; xyz < 3; ++xyz) std::cout<<std::setprecision(1)<<" "<<pfp.XYZ[startEnd][xyz] - truStart[xyz];
855  std::cout<<"\n";
856  }
857 */
858  hist.fPFPStartdX[truIndex]->Fill(pfp.XYZ[startEnd][0] - truStart[0]);
859  hist.fPFPStartdY[truIndex]->Fill(pfp.XYZ[startEnd][1] - truStart[1]);
860  hist.fPFPStartdZ[truIndex]->Fill(pfp.XYZ[startEnd][2] - truStart[2]);
861  Vector3_t truDir {{mcp->Px(), mcp->Py(), mcp->Pz()}};
862  SetMag(truDir, 1);
863  float dang = DeltaAngle(truDir, pfp.Dir[startEnd]);
864  while(dang > M_PI) dang -= twopi;
865  while(dang < -M_PI) dang += twopi;
866  hist.fPFPStartAngDiff[truIndex]->Fill(dang);
867  } // pfp
868 
869  StudyElectrons(hist);
870 
871  } // MatchTruth
872 
874  void TruthMatcher::MatchAndSum(const HistStuff& hist, const std::vector<unsigned int>& mcpSelect, const geo::TPCID& inTPCID)
875  {
876  // Match Tjs and PFParticles and accumulate performance statistics
877  if(mcpSelect.empty()) return;
878 
879  unsigned int tpc = inTPCID.TPC;
880  unsigned int cstat = inTPCID.Cryostat;
881 
882  // get the hits associated with all MCParticles in mcpSelect
883  std::vector<std::vector<unsigned int>> mcpHits(mcpSelect.size());
884  for(unsigned short isel = 0; isel < mcpSelect.size(); ++isel) {
885  unsigned int mcpIndex = mcpSelect[isel];
886  for(unsigned int iht = 0; iht < tjs.fHits.size(); ++iht) {
887  if(tjs.fHits[iht].MCPartListIndex != mcpIndex) continue;
888  if(tjs.fHits[iht].ArtPtr->WireID().TPC != tpc) continue;
889  if(tjs.fHits[iht].ArtPtr->WireID().Cryostat != cstat) continue;
890  mcpHits[isel].push_back(iht);
891  } // iht
892  } // mcpIndex
893 
894  // get the hits in all Tjs in this TPCID
895  std::vector<std::vector<unsigned int>> tjHits(tjs.allTraj.size());
896  for(unsigned short itj = 0; itj < tjs.allTraj.size(); ++itj) {
897  auto& tj = tjs.allTraj[itj];
898  if(tj.AlgMod[kKilled]) continue;
899  if(DecodeCTP(tj.CTP).TPC != tpc) continue;
900  tjHits[itj] = PutTrajHitsInVector(tj, kUsedHits);
901  tj.MCPartListIndex = UINT_MAX;
902  } // itj
903 
904  // match them up
905  for(unsigned short isel = 0; isel < mcpSelect.size(); ++isel) {
906  if(mcpHits[isel].empty()) continue;
907  unsigned int mcpIndex = mcpSelect[isel];
908  auto& mcp = tjs.MCPartList[mcpIndex];
909  unsigned short pdgIndex = PDGCodeIndex(tjs, mcp->PdgCode());
910  if(pdgIndex > 4) continue;
911  std::string particleName = "Other";
912  int pdg = abs(mcp->PdgCode());
913  if(pdg == 11) particleName = "Electron";
914  if(pdg == 22) particleName = "Photon";
915  if(pdg == 13) particleName = "Muon";
916  if(pdg == 211) particleName = "Pion";
917  if(pdg == 321) particleName = "Kaon";
918  if(pdg == 2212) particleName = "Proton";
919  if(particleName == "Other") particleName = "PDG_" + std::to_string(pdg);
920  float TMeV = 1000 * (mcp->E() - mcp->Mass());
921  for(unsigned short plane = 0; plane < tjs.NumPlanes; ++plane) {
922  // put the mcpHits (which are already in this TPCID) in this plane in a vector
923  std::vector<unsigned int> mcpPlnHits;
924  for(auto iht : mcpHits[isel]) {
925  auto& hit = tjs.fHits[iht];
926  if(hit.ArtPtr->WireID().Plane == plane) mcpPlnHits.push_back(iht);
927  } // iht
928  // require 2 truth-matched hits
929  if(mcpPlnHits.size() < 2) continue;
930  if((float)mcpPlnHits.size() >= tjs.MatchTruth[3]) ++nLongInPln;
931  TSums[pdgIndex] += TMeV;
932  ++EPCnts[pdgIndex];
933  CTP_t inCTP = EncodeCTP(cstat, tpc, plane);
934  unsigned short mtj = USHRT_MAX;
935  float maxEP = 0;
936  for(unsigned short itj = 0; itj < tjs.allTraj.size(); ++itj) {
937  // No hits in this TPC and plane
938  if(tjHits[itj].empty()) continue;
939  auto& tj = tjs.allTraj[itj];
940  // wrong CTP?
941  if(tj.CTP != inCTP) continue;
942  // make a list of hits that are common
943  auto shared = SetIntersection(mcpPlnHits, tjHits[itj]);
944  if(shared.empty()) continue;
945  float eff = (float)shared.size() / (float)mcpPlnHits.size();
946  float pur = (float)shared.size() / (float)tjHits[itj].size();
947  float ep = eff * pur;
948  // temp for checking
949  if(pdg != 11) {
950  hist.fEff->Fill(eff);
951  hist.fPur->Fill(pur);
952  }
953  // Replace a previously made poorer match with a better one?
954  if(tj.MCPartListIndex != UINT_MAX && ep > tj.EffPur) {
955  tj.EffPur = ep;
956  tj.MCPartListIndex = mcpIndex;
957  }
958  if(ep > maxEP) {
959  maxEP = ep;
960  mtj = itj;
961  }
962  } // itj
963  if(mtj == USHRT_MAX) {
964  // failed to match MCParticle to a trajectory
965  // Enter 0 in the profile histogram
966  hist.fEP_T[pdgIndex]->Fill(TMeV, 0);
967  if((float)mcpPlnHits.size() > tjs.MatchTruth[3]) {
968  ++nBadEP;
969  mf::LogVerbatim myprt("TC");
970  myprt<<particleName<<" BadEP TMeV "<<(int)TMeV<<" No matched trajectory to isel "<<isel;
971  myprt<<" nTrue hits "<<mcpPlnHits.size();
972  myprt<<" extent "<<PrintHit(tjs.fHits[mcpPlnHits[0]])<<"-"<<PrintHit(tjs.fHits[mcpPlnHits[mcpPlnHits.size() - 1]]);
973  myprt<<" events processed "<<tjs.EventsProcessed;
974  }
975  continue;
976  } // match failed
977  auto& tj = tjs.allTraj[mtj];
978  // don't clobber a better match
979  if(maxEP < tj.EffPur) continue;
980  tj.EffPur = maxEP;
981  tj.MCPartListIndex = mcpIndex;
982  EPTSums[pdgIndex] += TMeV * tj.EffPur;
983  hist.fEP_T[pdgIndex]->Fill(TMeV, tj.EffPur);
984  // ignore electrons
985  if(tj.EffPur < tjs.MatchTruth[2] && (float)mcpPlnHits.size() >= tjs.MatchTruth[3] && pdgIndex > 0) {
986  ++nBadEP;
987  mf::LogVerbatim myprt("TC");
988  myprt<<particleName<<" BadEP: "<<std::fixed<<std::setprecision(2)<<tj.EffPur;
989  myprt<<" mcpIndex "<<mcpIndex;
990  myprt<<" TMeV "<<(int)TMeV<<" MCP hits "<<mcpPlnHits.size();
991  myprt<<" extent "<<PrintHit(tjs.fHits[mcpPlnHits[0]])<<"-"<<PrintHit(tjs.fHits[mcpPlnHits[mcpPlnHits.size() - 1]]);
992  myprt<<" T"<<tj.ID;
993  myprt<<" Algs";
994  for(unsigned short ib = 0; ib < AlgBitNames.size(); ++ib) if(tj.AlgMod[ib]) myprt<<" "<<AlgBitNames[ib];
995  myprt<<" events processed "<<tjs.EventsProcessed;
996  } // print BadEP
997  // badep
998  } // plane
999  } // isel
1000 
1001  // Calculate PFParticle efficiency and purity
1002  std::vector<std::vector<unsigned int>> pfpHits;
1003  if(!tjs.pfps.empty()) {
1004  pfpHits.resize(tjs.pfps.size());
1005  // get the hits in all pfparticles in this TPCID
1006  for(unsigned short ipfp = 0; ipfp < tjs.pfps.size(); ++ipfp) {
1007  auto& pfp = tjs.pfps[ipfp];
1008  if(pfp.ID == 0) continue;
1009  // in the right TPCID?
1010  if(pfp.TPCID != inTPCID) continue;
1011  // ignore the neutrino PFParticle and true photons
1012  if(pfp.PDGCode == 14 || pfp.PDGCode == 12 || pfp.PDGCode == 22) continue;
1013  pfp.MCPartListIndex = UINT_MAX;
1014  for(auto& tjid : pfp.TjIDs) {
1015  unsigned short itj = tjid - 1;
1016  pfpHits[ipfp].insert(pfpHits[ipfp].end(), tjHits[itj].begin(), tjHits[itj].end());
1017  } // tj
1018  } // ipfp
1019  } // pfps exist
1020 
1021  // match them up
1022  for(unsigned short isel = 0; isel < mcpSelect.size(); ++isel) {
1023  unsigned short mpfp = USHRT_MAX;
1024  float maxEP = 0;
1025  unsigned int mcpIndex = mcpSelect[isel];
1026  if(mcpHits[isel].empty()) continue;
1027  auto& mcp = tjs.MCPartList[mcpIndex];
1028  float TMeV = 1000 * (mcp->E() - mcp->Mass());
1029  MCP_TSum += TMeV;
1030  // Performance reconstructing long muons, pions, kaons and protons
1031  unsigned short pdgIndex = PDGCodeIndex(tjs, mcp->PdgCode());
1032  bool longMCP = (pdgIndex > 0 && pdgIndex < 5 && (float)mcpHits[isel].size() >= 2 * tjs.MatchTruth[3]);
1033  if(longMCP) ++nLongMCP;
1034  for(unsigned short ipfp = 0; ipfp < tjs.pfps.size(); ++ipfp) {
1035  auto& pfp = tjs.pfps[ipfp];
1036  if(pfp.ID == 0) continue;
1037  // in the right TPCID?
1038  if(pfp.TPCID != inTPCID) continue;
1039  // not enough hits?
1040  if(pfpHits[ipfp].empty()) continue;
1041  auto shared = SetIntersection(mcpHits[isel], pfpHits[ipfp]);
1042  if(shared.empty()) continue;
1043  float eff = (float)shared.size() / (float)mcpHits[isel].size();
1044  float pur = (float)shared.size() / (float)pfpHits[ipfp].size();
1045  float ep = eff * pur;
1046  // Replace a previously made poorer match with a better one?
1047  if(pfp.MCPartListIndex != UINT_MAX && ep > pfp.EffPur) {
1048  pfp.EffPur = ep;
1049  pfp.MCPartListIndex = mcpIndex;
1050  }
1051  if(ep > maxEP) {
1052  maxEP = ep;
1053  mpfp = ipfp;
1054  }
1055  } // ipfp
1056  if(mpfp == USHRT_MAX) {
1057  if(TMeV > 30) {
1058  mf::LogVerbatim myprt("TC");
1059  myprt<<"BadPFP: MCParticle "<<mcpSelect[isel]<<" w PDGCode "<<mcp->PdgCode()<<" T "<<(int)TMeV<<" not reconstructed.";
1060  myprt<<" matched Tjs:";
1061  for(auto& tj : tjs.allTraj) {
1062  if(tj.AlgMod[kKilled]) continue;
1063  if(tj.MCPartListIndex == mcpIndex) myprt<<" "<<tj.ID<<" EP "<<std::fixed<<std::setprecision(2)<<tj.EffPur;
1064  } // tj
1065  myprt<<" events processed "<<tjs.EventsProcessed;
1066  } // TMeV > 30
1067  continue;
1068  }
1069  auto& pfp = tjs.pfps[mpfp];
1070  if(maxEP < pfp.EffPur) continue;
1071  pfp.EffPur = maxEP;
1072  pfp.MCPartListIndex = mcpIndex;
1073  MCP_EPTSum += TMeV * maxEP;
1074  ++MCP_PFP_Cnt;
1075  if(longMCP && maxEP > 0.8) ++nGoodLongMCP;
1076  } // isel
1077 
1078  MCP_Cnt += mcpSelect.size();
1079 
1080  } // MatchAndSum
1081 
1083  void TruthMatcher::PrintResults(int eventNum) const
1084  {
1085  // Print performance metrics for each selected event
1086  if(!tjs.SelectEvent) return;
1087 
1088  mf::LogVerbatim myprt("TC");
1089  myprt<<"Evt "<<eventNum;
1090  float sum = 0;
1091  float sumt = 0;
1092  for(unsigned short pdgIndex = 0; pdgIndex < TSums.size(); ++pdgIndex) {
1093  if(TSums[pdgIndex] == 0) continue;
1094  if(pdgIndex == 0) myprt<<" El";
1095  if(pdgIndex == 1) myprt<<" Mu";
1096  if(pdgIndex == 2) myprt<<" Pi";
1097  if(pdgIndex == 3) myprt<<" K";
1098  if(pdgIndex == 4) myprt<<" P";
1099  float ave = EPTSums[pdgIndex] / (float)TSums[pdgIndex];
1100  myprt<<" "<<std::fixed<<std::setprecision(2)<<ave;
1101 // myprt<<" "<<EPCnts[pdgIndex];
1102  if(pdgIndex > 0) {
1103  sum += TSums[pdgIndex];
1104  sumt += EPTSums[pdgIndex];
1105  }
1106  } // pdgIndex
1107  if(sum > 0) myprt<<" MuPiKP "<<std::fixed<<std::setprecision(2)<<sumt / sum;
1108  myprt<<" BadEP "<<nBadEP;
1109  if(nLongInPln > 0) {
1110  float longGood = 1 - (float)nBadEP / (float)nLongInPln;
1111  myprt<<" longGood "<<std::fixed<<std::setprecision(2)<<longGood;
1112  }
1113  if(MCP_TSum > 0) {
1114  // PFParticle statistics
1115  float ep = MCP_EPTSum / MCP_TSum;
1116  myprt<<" MCP cnt "<<(int)MCP_Cnt<<" PFP "<<std::fixed<<std::setprecision(2)<<ep;
1117  }
1118  if(Prim_TSum > 0) {
1119  float ep = Prim_EPTSum / Prim_TSum;
1120  myprt<<" PrimPFP "<<std::fixed<<std::setprecision(2)<<ep;
1121  }
1122  if(nLongMCP > 0) {
1123  float longGood = (float)nGoodLongMCP / (float)nLongMCP;
1124  myprt<<" longGood "<<std::fixed<<std::setprecision(2)<<longGood;
1125  }
1126  if(TruVxCounts[1] > 0) {
1127  // True vertex is reconstructable
1128  float frac = (float)TruVxCounts[2] / (float)TruVxCounts[1];
1129  myprt<<" NuVx correct "<<std::fixed<<std::setprecision(2)<<frac;
1130  }
1131 
1132  } // PrintResults
1133 
1135  bool TruthMatcher::CanReconstruct(unsigned int mcpIndex, unsigned short nDimensions, const geo::TPCID& inTPCID)
1136  {
1137  // returns true if the MCParticle can be reconstructed in nDimensions
1138  if(mcpIndex > tjs.MCPartList.size() - 1) return false;
1139  if(nDimensions < 2 || nDimensions > 3) return false;
1140 
1141  std::vector<unsigned short> cntInPln(tjs.NumPlanes);
1142  for(auto& hit : tjs.fHits) {
1143  if(hit.ArtPtr->WireID().TPC != inTPCID.TPC) continue;
1144  if(hit.ArtPtr->WireID().Cryostat != inTPCID.Cryostat) continue;
1145  if(hit.MCPartListIndex == mcpIndex) ++cntInPln[hit.ArtPtr->WireID().Plane];
1146  } // hit
1147  unsigned short nPlnOK = 0;
1148  // Require at least 2 truth-matched hits in a plane
1149  for(unsigned short plane = 0; plane < tjs.NumPlanes; ++plane) if(cntInPln[plane] > 1) ++nPlnOK;
1150  if(nPlnOK < nDimensions - 1) return false;
1151  return true;
1152  } // CanReconstruct
1153 
1154 
1156  std::vector<unsigned int> TruthMatcher::PutMCPHitsInVector(unsigned int mcpIndex, CTP_t inCTP)
1157  {
1158  // put the hits matched to the MCParticle into a vector in the requested CTP
1159  std::vector<unsigned int> hitVec;
1160  if(mcpIndex > tjs.MCPartList.size() - 1) return hitVec;
1161  geo::PlaneID planeID = DecodeCTP(inCTP);
1162  for(unsigned int iht = 0; iht < tjs.fHits.size(); ++iht) {
1163  auto& hit = tjs.fHits[iht];
1164  if(hit.MCPartListIndex != mcpIndex) continue;
1165  if(hit.ArtPtr->WireID().Plane != planeID.Plane) continue;
1166  if(hit.ArtPtr->WireID().TPC != planeID.TPC) continue;
1167  if(hit.ArtPtr->WireID().Cryostat != planeID.Cryostat) continue;
1168  hitVec.push_back(iht);
1169  } // iht
1170  return hitVec;
1171  } // PutMCPHitsInVector
1172 
1174  void MCParticleListUtils::MakeTruTrajPoint(unsigned int MCParticleListIndex, TrajPoint& tp)
1175  {
1176  // Creates a trajectory point at the start of the MCParticle with index MCParticleListIndex. The
1177  // projected length of the MCParticle in the plane coordinate system is stored in TruTp.Delta.
1178  // The calling function should specify the CTP in which this TP resides.
1179 
1180  if(MCParticleListIndex > tjs.MCPartList.size() - 1) return;
1181 
1182  const simb::MCParticle* mcp = tjs.MCPartList[MCParticleListIndex];
1183 
1184  Point3_t pos {{mcp->Vx(), mcp->Vy(), mcp->Vz()}};
1185  Vector3_t dir {{mcp->Px(), mcp->Py(), mcp->Pz()}};
1186  SetMag(dir, 1);
1187  tp = MakeBareTP(tjs, pos, dir, tp.CTP);
1188  } // MakeTruTrajPoint
1189 
1191  ShowerStruct3D MCParticleListUtils::MakeCheatShower(unsigned int mcpIndex, Point3_t primVx, int& truParentPFP)
1192  {
1193  // Make a 3D shower for an electron or photon MCParticle in the current TPC.
1194  // The shower ID is set to 0 if there is a failure. The ID of the most likely parent pfp is returned
1195  // as well - the one that is closest to the primary vertex position
1196  auto ss3 = CreateSS3(tjs, tjs.TPCID);
1197  truParentPFP = 0;
1198  // save the ID
1199  int goodID = ss3.ID;
1200  // set it to the failure state
1201  ss3.ID = 0;
1202  ss3.Cheat = true;
1203  if(mcpIndex > tjs.MCPartList.size() - 1) return ss3;
1204  auto& mcp = tjs.MCPartList[mcpIndex];
1205  int pdg = abs(mcp->PdgCode());
1206  if(!(pdg == 11 || pdg == 111)) return ss3;
1208  int eveID = mcp->TrackId();
1209  float showerEnergy = 1000 * mcp->E();
1210  float shMaxAlong = 7.0 * log(showerEnergy / 15);
1211  std::cout<<"MCS: showerEnergy "<<(int)showerEnergy<<" MeV. shMaxAlong "<<shMaxAlong<<" cm\n";
1212 
1213  // put the shower hits in two vectors. One for the InShower hits
1214  // and one for the shower parent
1215  std::vector<std::vector<unsigned int>> showerHits(tjs.NumPlanes);
1216  std::vector<std::vector<unsigned int>> parentHits(tjs.NumPlanes);
1217  unsigned short nsh = 0;
1218  unsigned short npar = 0;
1219  unsigned int cstat = tjs.TPCID.Cryostat;
1220  unsigned int tpc = tjs.TPCID.TPC;
1221  for(unsigned int iht = 0; iht < tjs.fHits.size(); ++iht) {
1222  auto& hit = tjs.fHits[iht];
1223  if(hit.MCPartListIndex > tjs.MCPartList.size()-1) continue;
1224  if(hit.ArtPtr->WireID().TPC != tpc) continue;
1225  if(hit.ArtPtr->WireID().Cryostat != cstat) continue;
1226  if(hit.Integral <= 0) continue;
1227  // require that it be used in a tj to ignore hits that are far
1228  // from the shower
1229  if(hit.InTraj <= 0) continue;
1230  auto& shmcp = tjs.MCPartList[hit.MCPartListIndex];
1231  // look for an electron
1232  if(abs(shmcp->PdgCode()) != 11) continue;
1233  // with eve ID = the primary electron being considered
1234  int eid = pi_serv->ParticleList().EveId(shmcp->TrackId());
1235  if(eid != eveID) continue;
1236  if(shmcp->TrackId() == eid) {
1237  // store the shower parent hit
1238  parentHits[hit.ArtPtr->WireID().Plane].push_back(iht);
1239  ++npar;
1240  } else {
1241  // store the InShower hit
1242  showerHits[hit.ArtPtr->WireID().Plane].push_back(iht);
1243  ++nsh;
1244  }
1245  } // iht
1246  // require more hits in the shower than in the parent
1247  if(npar > nsh) return ss3;
1248 /*
1249  for(unsigned short plane = 0; plane < tjs.NumPlanes; ++plane) {
1250  std::cout<<" plane "<<plane<<" shower hits "<<showerHits[plane].size();
1251  std::cout<<" parentHits "<<parentHits[plane].size()<<"\n";
1252  } // plane
1253 */
1254  // create 2D cheat showers in each plane
1255  std::vector<int> dummy_tjlist;
1256  for(unsigned short plane = 0; plane < tjs.NumPlanes; ++plane) {
1257  CTP_t inCTP = EncodeCTP(cstat, tpc, plane);
1258  auto ss = CreateSS(tjs, inCTP, dummy_tjlist);
1259  // fill the ShPts
1260  ss.ShPts.resize(showerHits[plane].size());
1261  for(unsigned short iht = 0; iht < showerHits[plane].size(); ++iht) {
1262  auto& hit = tjs.fHits[showerHits[plane][iht]];
1263  ss.ShPts[iht].HitIndex = showerHits[plane][iht];
1264  ss.ShPts[iht].TID = 0;
1265  ss.ShPts[iht].Chg = hit.Integral;
1266  ss.ShPts[iht].Pos[0] = hit.ArtPtr->WireID().Wire;
1267  ss.ShPts[iht].Pos[1] = hit.PeakTime * tjs.UnitsPerTick;
1268  } // iht
1269  ss.SS3ID = goodID;
1270  if(!UpdateShower("MCS", tjs, ss, true)) {
1271  std::cout<<"Failed to update 2S"<<ss.ID<<"\n";
1272  return ss3;
1273  } // UpdateShower failed
1274  // remove shower points that are far away
1275  std::vector<ShowerPoint> spts;
1276  for(auto& spt : ss.ShPts) {
1277  // don't make a tight cut right now. The shower parameterization isn't great
1278  double along = tjs.WirePitch * spt.RotPos[0];
1279  float tau = along / shMaxAlong;
1280 // auto& hit = tjs.fHits[spt.HitIndex];
1281 // if(ss.ID == 2) std::cout<<"chk "<<PrintHit(hit)<<" along "<<(int)along<<" tau "<<std::fixed<<std::setprecision(1)<<tau<<"\n";
1282  if(tau > -1 && tau < 2) {
1283  spts.push_back(spt);
1284  } else {
1285 // std::cout<<"Skip "<<PrintHit(hit)<<" along "<<(int)along<<"\n";
1286  }
1287  } // spt
1288  std::cout<<" 2S"<<ss.ID<<" shpts "<<ss.ShPts.size()<<" new "<<spts.size()<<"\n";
1289  ss.ShPts = spts;
1290  ss.NeedsUpdate = true;
1291  if(!UpdateShower("MCS", tjs, ss, true)) {
1292  std::cout<<"Failed to update 2S"<<ss.ID<<"\n";
1293  return ss3;
1294  } // UpdateShower failed
1295  if(!StoreShower("MCS", tjs, ss)) {
1296  std::cout<<"Failed to store 2S"<<ss.ID<<"\n";
1297  return ss3;
1298  } // UpdateShower failed
1299  ss3.CotIDs.push_back(ss.ID);
1300  } // plane
1301  ss3.ID = goodID;
1302  if(!UpdateShower("MCS", tjs, ss3, true)) {
1303  std::cout<<"SS3 Failed...\n";
1304  ss3.ID = 0;
1305  return ss3;
1306  }
1307  // success
1308 
1309  // now look for a parent pfp
1310  std::vector<std::pair<int, int>> pcnt;
1311  for(unsigned short plane = 0; plane < tjs.NumPlanes; ++plane) {
1312  for(auto iht : parentHits[plane]) {
1313  auto& hit = tjs.fHits[iht];
1314  if(hit.InTraj <= 0) continue;
1315  auto& tj = tjs.allTraj[hit.InTraj - 1];
1316  if(!tj.AlgMod[kMat3D]) continue;
1317  auto TInP = GetAssns(tjs, "T", tj.ID, "P");
1318  if(TInP.empty()) continue;
1319  auto& pfp = tjs.pfps[TInP[0] - 1];
1320  unsigned short indx = 0;
1321  for(indx = 0; indx < pcnt.size(); ++indx) if(pfp.ID == pcnt[indx].first) break;
1322  if(indx == pcnt.size()) {
1323  pcnt.push_back(std::make_pair(pfp.ID, 1));
1324  } else {
1325  ++pcnt[indx].second;
1326  }
1327  } // iht
1328  } // plane
1329  float close = 5;
1330  for(auto pcn : pcnt) {
1331  if(pcn.second < 6) continue;
1332  auto& pfp = tjs.pfps[pcn.first - 1];
1333  for(unsigned short end = 0; end < 2; ++end) {
1334  float sep = PosSep(pfp.XYZ[end], primVx);
1335  if(sep > close) continue;
1336  close = sep;
1337  truParentPFP = pfp.ID;
1338  }
1339  } // pcn
1340  std::cout<<"truParent P"<<truParentPFP<<" close "<<std::fixed<<std::setprecision(2)<<close<<"\n";
1341 
1342  return ss3;
1343  } // MakeCheatShower
1344 
1347  {
1348  // returns the SCE corrected start position of a primary electron
1349  if(tjs.MCPartList.empty()) return false;
1351  for(unsigned int part = 0; part < tjs.MCPartList.size(); ++part) {
1352  auto& mcp = tjs.MCPartList[part];
1353  // require electron
1354  if(abs(mcp->PdgCode()) != 11) continue;
1355  int eveID = pi_serv->ParticleList().EveId(mcp->TrackId());
1356  if(mcp->TrackId() != eveID) continue;
1357  start = {{mcp->Vx(), mcp->Vy(), mcp->Vz()}};
1358  auto const* SCE = lar::providerFrom<spacecharge::SpaceChargeService>();
1359  geo::Vector_t posOffsets = SCE->GetPosOffsets({start[0], start[1], start[2]});
1360  posOffsets.SetX(-posOffsets.X());
1361  start[0] += posOffsets.X();
1362  start[1] += posOffsets.Y();
1363  start[2] += posOffsets.Z();
1364  dir = {{mcp->Px(), mcp->Py(), mcp->Pz()}};
1365  SetMag(dir, 1);
1366  energy = 1000 * (mcp->E() - mcp->Mass());
1367  return true;
1368  } // part
1369  return false;
1370  } // PrimaryElectronStart
1371 
1372 
1375  {
1376  // Returns the ID of a pfp that has hits whose eve ID is a primary electron
1377  // and is the closest (< 5 WSE units) to the primary electron start
1378  if(tjs.MCPartList.empty()) return 0;
1380  TruthMatcher tm{tjs};
1381  for(unsigned int part = 0; part < tjs.MCPartList.size(); ++part) {
1382  auto& mcp = tjs.MCPartList[part];
1383  // require electron
1384  if(abs(mcp->PdgCode()) != 11) continue;
1385  int eveID = pi_serv->ParticleList().EveId(mcp->TrackId());
1386  if(mcp->TrackId() != eveID) continue;
1387  Point3_t start = {{mcp->Vx(), mcp->Vy(), mcp->Vz()}};
1388  auto const* SCE = lar::providerFrom<spacecharge::SpaceChargeService>();
1389  geo::Vector_t posOffsets = SCE->GetPosOffsets({start[0], start[1], start[2]});
1390  posOffsets.SetX(-posOffsets.X());
1391  start[0] += posOffsets.X();
1392  start[1] += posOffsets.Y();
1393  start[2] += posOffsets.Z();
1394  // make a list of pfps that use tjs that use these hits
1395  std::vector<int> pfplist;
1396  for(unsigned short plane = 0; plane < tjs.NumPlanes; ++plane) {
1397  CTP_t inCTP = EncodeCTP(tpcid.Cryostat, tpcid.TPC, plane);
1398  std::vector<unsigned int> mcphits = tm.PutMCPHitsInVector(part, inCTP);
1399  if(mcphits.empty()) continue;
1400  for(auto iht : mcphits) {
1401  auto& hit = tjs.fHits[iht];
1402  if(hit.InTraj <= 0) continue;
1403  // require that the tj is 3D-matched
1404  auto& tj = tjs.allTraj[hit.InTraj - 1];
1405  if(!tj.AlgMod[kMat3D]) continue;
1406  // find the number of hits that are in mcphits
1407  auto tjhits = PutTrajHitsInVector(tj, kUsedHits);
1408  auto shared = SetIntersection(tjhits, mcphits);
1409  if(shared.size() < 10) continue;
1410  // find out what pfp it is used in.
1411  auto TInP = GetAssns(tjs, "T", tj.ID, "P");
1412  if(TInP.size() != 1) continue;
1413  int pid = TInP[0];
1414  if(std::find(pfplist.begin(), pfplist.end(), pid) == pfplist.end()) pfplist.push_back(pid);
1415  } // iht
1416  } // plane
1417  if(pfplist.empty()) return 0;
1418  // Use the one that is closest to the true start position, not the
1419  // one that has the most matching hits. Electrons are likely to be
1420  // poorly reconstructed
1421  int pfpid = 0;
1422  float close = 1E6;
1423  for(auto pid : pfplist) {
1424  auto& pfp = tjs.pfps[pid - 1];
1425  unsigned short nearEnd = 1 - FarEnd(tjs, pfp, start);
1426  float sep = PosSep2(pfp.XYZ[nearEnd], start);
1427  if(sep > close) continue;
1428  close = sep;
1429  pfpid = pid;
1430  }
1431  return pfpid;
1432  } // part
1433  return 0;
1434  } // PrimaryElectronPFPID
1435 
1438  {
1439  // returns the ID of a tj in inCTP that is closest to the start of a primary electron
1440  if(tjs.MCPartList.empty()) return 0;
1442  for(unsigned int part = 0; part < tjs.MCPartList.size(); ++part) {
1443  auto& mcp = tjs.MCPartList[part];
1444  // require electron
1445  if(abs(mcp->PdgCode()) != 11) continue;
1446  int eveID = pi_serv->ParticleList().EveId(mcp->TrackId());
1447  if(mcp->TrackId() != eveID) continue;
1448  return MCParticleStartTjID(part, inCTP);
1449  } // part
1450  return 0;
1451  } // PrimaryElectronTjID
1452 
1454  int MCParticleListUtils::MCParticleStartTjID(unsigned int mcpIndex, CTP_t inCTP)
1455  {
1456  // Finds the trajectory that has hits matched to the MC Particle and is the closest to the
1457  // MCParticle start vertex
1458 
1459  if(mcpIndex > tjs.MCPartList.size() - 1) return 0;
1460 
1461  geo::PlaneID planeID = DecodeCTP(inCTP);
1462 
1463  // tj ID and occurrence count
1464  std::vector<std::array<int, 2>> t_cnt;
1465  for(auto hit : tjs.fHits) {
1466  if(hit.InTraj <= 0) continue;
1467  if(hit.MCPartListIndex != mcpIndex) continue;
1468  if(hit.ArtPtr->WireID().TPC != planeID.TPC) continue;
1469  if(hit.ArtPtr->WireID().Plane != planeID.Plane) continue;
1470  if(hit.ArtPtr->WireID().Cryostat != planeID.Cryostat) continue;
1471  unsigned short indx = 0;
1472  for(indx = 0; indx < t_cnt.size(); ++indx) if(t_cnt[indx][0] == hit.InTraj) break;
1473  if(indx == t_cnt.size()) {
1474  // didn't find the traj ID in t_cnt so add it
1475  std::array<int, 2> tmp;
1476  tmp[0] = hit.InTraj;
1477  tmp[1] = 1;
1478  t_cnt.push_back(tmp);
1479  } else {
1480  // count occurrences of this tj ID
1481  ++t_cnt[indx][1];
1482  }
1483  } // hit
1484 
1485  if(t_cnt.empty()) return 0;
1486  unsigned short occMax = 0;
1487  unsigned short occIndx = 0;
1488  for(unsigned short ii = 0; ii < t_cnt.size(); ++ii) {
1489  auto& tcnt = t_cnt[ii];
1490  if(tcnt[1] < occMax) continue;
1491  occMax = tcnt[1];
1492  occIndx = ii;
1493  } // tcnt
1494  return t_cnt[occIndx][0];
1495 
1496  } // MCParticleStartTj
1497 
1500  {
1501  // Returns the MCParticle index that best matches the hits used in the Tp
1502  if(tjs.MCPartList.empty()) return UINT_MAX;
1503  if(tp.Chg <= 0) return UINT_MAX;
1504  std::vector<unsigned int> mcpIndex;
1505  std::vector<unsigned short> mcpCnt;
1506  for(unsigned short ii = 0; ii < tp.Hits.size(); ++ii) {
1507  if(!tp.UseHit[ii]) continue;
1508  unsigned int mcpi = tjs.fHits[tp.Hits[ii]].MCPartListIndex;
1509  if(mcpi == UINT_MAX) continue;
1510  unsigned short indx = 0;
1511  for(indx = 0; indx < mcpIndex.size(); ++indx) if(mcpi == mcpIndex[indx]) break;
1512  if(indx == mcpIndex.size()) {
1513  // not in the list so add it
1514  mcpIndex.push_back(mcpi);
1515  mcpCnt.push_back(1);
1516  } else {
1517  ++mcpCnt[indx];
1518  }
1519  } // ii
1520  if(mcpIndex.empty()) return UINT_MAX;
1521  if(mcpIndex.size() == 1) return mcpIndex[0];
1522  unsigned int indx = 0;
1523  unsigned short maxCnt = 0;
1524  for(unsigned short ii = 0; ii < mcpIndex.size(); ++ii) {
1525  if(mcpCnt[ii] > maxCnt) {
1526  maxCnt = mcpCnt[ii];
1527  indx = mcpIndex[ii];
1528  }
1529  } // ii
1530  return indx;
1531  } // GetMCPartListIndex
1532 
1534  unsigned int MCParticleListUtils::GetMCPartListIndex(const ShowerStruct& ss, unsigned short& nTruHits)
1535  {
1536  // Returns the index of the MCParticle that has the most number of matches
1537  // to the hits in this shower
1538 
1539  if(tjs.MCPartList.empty()) return UINT_MAX;
1540  if(ss.TjIDs.empty()) return UINT_MAX;
1541 
1542  std::vector<unsigned int> pListCnt(tjs.MCPartList.size());
1543 
1544  for(auto& tjid : ss.TjIDs) {
1545  Trajectory& tj = tjs.allTraj[tjid - 1];
1546  for(auto& tp : tj.Pts) {
1547  for(unsigned short ii = 0; ii < tp.Hits.size(); ++ii) {
1548  if(!tp.UseHit[ii]) continue;
1549  unsigned int iht = tp.Hits[ii];
1550  // ignore unmatched hits
1551  if(tjs.fHits[iht].MCPartListIndex > tjs.MCPartList.size() - 1) continue;
1552  ++pListCnt[tjs.fHits[iht].MCPartListIndex];
1553  } // ii
1554  } // pt
1555  } // tjid
1556 
1557  unsigned int pIndex = UINT_MAX;
1558  nTruHits = 0;
1559  for(unsigned short ii = 0; ii < pListCnt.size(); ++ii) {
1560  if(pListCnt[ii] > nTruHits) {
1561  nTruHits = pListCnt[ii];
1562  pIndex = ii;
1563  }
1564  } // ii
1565 
1566  return pIndex;
1567 
1568  } // GetMCPartListIndex
1569 
1571  unsigned int MCParticleListUtils::GetMCPartListIndex(const Trajectory& tj, unsigned short& nTruHits)
1572  {
1573  // Returns the index of the MCParticle that has the most number of matches
1574  // to the hits in this trajectory
1575 
1576  if(tjs.MCPartList.empty()) return UINT_MAX;
1577 
1578  // Check all hits associated with this Tj
1579  std::vector<unsigned int> pListCnt(tjs.MCPartList.size());
1580 
1581  for(auto& tp : tj.Pts) {
1582  for(unsigned short ii = 0; ii < tp.Hits.size(); ++ii) {
1583  if(!tp.UseHit[ii]) continue;
1584  unsigned int iht = tp.Hits[ii];
1585  // ignore unmatched hits
1586  if(tjs.fHits[iht].MCPartListIndex > tjs.MCPartList.size() - 1) continue;
1587  ++pListCnt[tjs.fHits[iht].MCPartListIndex];
1588  } // ii
1589  } // pt
1590 
1591  unsigned int pIndex = UINT_MAX;
1592  nTruHits = 0;
1593  for(unsigned short ii = 0; ii < pListCnt.size(); ++ii) {
1594  if(pListCnt[ii] > nTruHits) {
1595  nTruHits = pListCnt[ii];
1596  pIndex = ii;
1597  }
1598  } // ii
1599 
1600  return pIndex;
1601 
1602  } // GetMCPartListIndex
1603 
1604 } // namespace tca
double E(const int i=0) const
Definition: MCParticle.h:237
std::bitset< 16 > UseHit
Definition: DataStructs.h:163
float fShEnergy
Definition: TCHist.h:86
float MCP_TSum
Definition: TCTruth.h:61
TH1F * fTruT[5]
Definition: TCHist.h:40
ROOT::Math::DisplacementVector3D< ROOT::Math::Cartesian3D< double >, ROOT::Math::GlobalCoordinateSystemTag > Vector_t
Type for representation of momenta in 3D space.
Definition: geo_vectors.h:167
MaybeLogger_< ELseverityLevel::ELsev_info, true > LogVerbatim
bool SelectEvent
number of points to find AveChg
Definition: DataStructs.h:533
int PdgCode() const
Definition: MCParticle.h:216
bool empty() const
Definition: ParticleList.h:314
float fPfpLen
Definition: TCHist.h:86
unsigned short nLongInPln
Definition: TCTruth.h:70
TProfile * fNearTj[5]
Definition: TCHist.h:70
double Py(const int i=0) const
Definition: MCParticle.h:235
const std::vector< std::string > AlgBitNames
Definition: DataStructs.cxx:4
CTP_t CTP
Cryostat, TPC, Plane code.
Definition: DataStructs.h:145
float fTrans
Definition: TCHist.h:86
double InShowerProbLong(double showerEnergy, double along)
Definition: TCShower.cxx:2118
std::array< double, 3 > Point3_t
Definition: DataStructs.h:35
TH1F * fPur
Definition: TCHist.h:81
const std::vector< sim::TrackIDE > HitToTrackIDEs(recob::Hit const &hit)
bool PrimaryElectronStart(Point3_t &start, Vector3_t &dir, float &energy)
Definition: TCTruth.cxx:1346
float fChgFrac
Definition: TCHist.h:86
TH1F * fElectronLike[5]
Definition: TCHist.h:28
unsigned int Nplanes() const
Number of planes in this tpc.
Definition: TPCGeo.h:145
std::vector< unsigned int > PutTrajHitsInVector(Trajectory const &tj, HitStatus_t hitRequest)
Definition: Utils.cxx:2283
int Mother() const
Definition: MCParticle.h:217
float ChgFracBetween(TjStuff &tjs, Point3_t pos1, Point3_t pos2, geo::TPCID tpcid)
Definition: PFPUtils.cxx:2860
unsigned short PDGCodeIndex(TjStuff &tjs, int PDGCode)
Definition: Utils.cxx:1783
void FindAlongTrans(Point3_t pos1, Vector3_t dir1, Point3_t pos2, Point2_t &alongTrans)
Definition: PFPUtils.cxx:2751
TH1F * fPFPStartEnd
Definition: TCHist.h:74
list_type::const_iterator const_iterator
Definition: ParticleList.h:132
simb::Origin_t Origin() const
Definition: MCTruth.h:71
double Mass() const
Definition: MCParticle.h:243
The data type to uniquely identify a Plane.
Definition: geo_types.h:250
Geometry information for a single TPC.
Definition: TPCGeo.h:37
std::vector< PFPStruct > pfps
Definition: DataStructs.h:505
enum simb::_ev_origin Origin_t
event origin types
std::vector< unsigned int > Hits
Definition: DataStructs.h:162
TTree * fShowerParentBkg
Definition: TCHist.h:84
double Px(const int i=0) const
Definition: MCParticle.h:234
std::vector< int > TjIDs
Definition: DataStructs.h:275
TProfile * fEP_T[5]
Definition: TCHist.h:67
TH2F * PDGCode_reco_true
Definition: TCHist.h:73
Float_t ss
Definition: plot.C:23
CryostatID_t Cryostat
Index of cryostat.
Definition: geo_types.h:130
int PrimaryElectronPFPID(const geo::TPCID &tpcid)
Definition: TCTruth.cxx:1374
std::vector< ShowerStruct3D > showers
Definition: DataStructs.h:508
Float_t tmp
Definition: plot.C:37
ShowerStruct3D MakeCheatShower(unsigned int mcpIndex, Point3_t primVx, int &truParentPFP)
Definition: TCTruth.cxx:1191
TH1F * fMomAsym[5]
Definition: TCHist.h:27
double DeltaAngle(const Vector3_t v1, const Vector3_t v2)
Definition: PFPUtils.cxx:1607
double ShowerEnergy(const ShowerStruct3D &ss3)
Definition: TCShower.cxx:4434
SigType_t SignalType(geo::PlaneID const &pid) const
Returns the type of signal on the channels of specified TPC plane.
int EveId(const int trackID) const
TH1F * fPFPStartdZ[5]
Definition: TCHist.h:77
std::string Process() const
Definition: MCParticle.h:219
unsigned short NumPlanes
Definition: DataStructs.h:475
int PrimaryElectronTjID(CTP_t inCTP)
Definition: TCTruth.cxx:1437
void Initialize()
Definition: TCTruth.cxx:20
IteratorBox< TPC_id_iterator,&GeometryCore::begin_TPC_id,&GeometryCore::end_TPC_id > IterateTPCIDs() const
Enables ranged-for loops on all TPC IDs of the detector.
void StudyShowerParents(HistStuff &hist)
Definition: TCTruth.cxx:235
std::vector< T > SetIntersection(const std::vector< T > &set1, const std::vector< T > &set2)
Definition: Utils.h:214
int MCParticleStartTjID(unsigned int MCParticleListIndex, CTP_t inCTP)
Definition: TCTruth.cxx:1454
unsigned short FarEnd(const TjStuff &tjs, const PFPStruct &pfp, const Point3_t &pos)
Definition: PFPUtils.cxx:2936
bool CanReconstruct(unsigned int mcpIndex, unsigned short nDimensions, const geo::TPCID &tpcid)
Definition: TCTruth.cxx:1135
int TrackId() const
Definition: MCParticle.h:214
TH2F * fElectronLike_Len[5]
Definition: TCHist.h:29
float Prim_EPTSum
Definition: TCTruth.h:66
TH1F * fPFPStartAngDiff[5]
Definition: TCHist.h:78
TjStuff & tjs
Definition: TCTruth.h:55
float MCP_EPTSum
Definition: TCTruth.h:62
std::vector< ShowerStruct > cots
Definition: DataStructs.h:506
float fAlong
Definition: TCHist.h:86
TH1F * fChgToMeV[3]
Definition: TCHist.h:31
unsigned int EventsProcessed
Definition: DataStructs.h:521
bool UpdateShower(std::string inFcnLabel, TjStuff &tjs, ShowerStruct &ss, bool prt)
Definition: TCShower.cxx:1104
std::array< float, 2 > Point2_t
Definition: DataStructs.h:37
const geo::GeometryCore * geom
Definition: DataStructs.h:526
Vector3_t PointDirection(const Point3_t p1, const Point3_t p2)
Definition: PFPUtils.cxx:1614
void MatchTrueHits()
Definition: TCTruth.cxx:31
single particles thrown at the detector
Definition: MCTruth.h:24
std::vector< TrajPoint > Pts
Trajectory points.
Definition: DataStructs.h:169
std::array< float, 5 > TSums
Definition: TCTruth.h:58
float fSep
Definition: TCHist.h:86
TrajPoint MakeBareTP(TjStuff &tjs, Point3_t &pos, Vector3_t &dir, CTP_t inCTP)
Definition: Utils.cxx:3435
std::array< float, 5 > EPTSums
Definition: TCTruth.h:59
std::vector< evd::details::RawDigitInfo_t >::const_iterator begin(RawDigitCacheDataClass const &cache)
TString part[npart]
Definition: Style.C:32
std::array< short, 5 > EPCnts
Definition: TCTruth.h:57
std::vector< Trajectory > allTraj
vector of all trajectories in each plane
Definition: DataStructs.h:493
double energy
Definition: plottest35.C:25
std::vector< unsigned int > PutMCPHitsInVector(unsigned int mcpIndex, CTP_t inCTP)
Definition: TCTruth.cxx:1156
iterator begin()
Definition: ParticleList.h:305
TH1F * fEff
Definition: TCHist.h:80
double PosSep2(const Point3_t &pos1, const Point3_t &pos2)
Definition: PFPUtils.cxx:1631
std::string PrintHit(const TCHit &hit)
Definition: Utils.cxx:4732
bool StoreShower(std::string inFcnLabel, TjStuff &tjs, ShowerStruct3D &ss3)
Definition: TCShower.cxx:4468
float fPfpEnergy
Definition: TCHist.h:86
unsigned short nLongMCP
Definition: TCTruth.h:71
float UnitsPerTick
scale factor from Tick to WSE equivalent units
Definition: DataStructs.h:468
The data type to uniquely identify a TPC.
Definition: geo_types.h:195
PlaneID_t Plane
Index of the plane within its TPC.
Definition: geo_types.h:258
TTree * fShowerParentSig
Definition: TCHist.h:83
float fInShwrProb
Definition: TCHist.h:86
void MatchTruth(const HistStuff &hist, bool fStudyMode)
Definition: TCTruth.cxx:459
unsigned short nBadEP
Definition: TCTruth.h:69
unsigned short nGoodLongMCP
Definition: TCTruth.h:72
bool SetMag(Vector3_t &v1, double mag)
Definition: PFPUtils.cxx:1641
View_t View(geo::PlaneID const &pid) const
Returns the view (wire orientation) on the channels of specified TPC plane.
geo::TPCID TPCID
Definition: DataStructs.h:469
void MatchAndSum(const HistStuff &hist, const std::vector< unsigned int > &mcpSelect, const geo::TPCID &inTPCID)
Definition: TCTruth.cxx:874
std::vector< float > MatchTruth
Match to MC truth.
Definition: DataStructs.h:518
Detector simulation of raw signals on wires.
short MCSMom(const TjStuff &tjs, const std::vector< int > &tjIDs)
Definition: Utils.cxx:2837
std::vector< Vtx3Store > vtx3
3D vertices
Definition: DataStructs.h:503
float fDang2
Definition: TCHist.h:86
double Vx(const int i=0) const
Definition: MCParticle.h:225
ShowerStruct3D CreateSS3(TjStuff &tjs, const geo::TPCID &tpcid)
Definition: TCShower.cxx:4553
raw::ChannelID_t PlaneWireToChannel(WireID const &wireid) const
Returns the ID of the TPC channel connected to the specified wire.
TH1F * fChgRMS[5]
Definition: TCHist.h:26
double DotProd(const Vector3_t &v1, const Vector3_t &v2)
Definition: PFPUtils.h:60
unsigned int CTP_t
Definition: DataStructs.h:41
std::string to_string(Flag_t< Storage > const flag)
Convert a flag into a stream (shows its index).
Definition: BitMask.h:187
std::vector< TCHit > fHits
Definition: DataStructs.h:495
in close()
TH2F * hist
Definition: plot.C:136
std::string PrintHitShort(const TCHit &hit)
Definition: Utils.cxx:4726
TDirectory * dir
Definition: macro.C:5
TProfile * fChgToMeV_Etru
Definition: TCHist.h:32
TH1F * fPFPStartdY[5]
Definition: TCHist.h:76
float WirePitch
Definition: DataStructs.h:482
TH1F * fPFPStartdX[5]
Definition: TCHist.h:75
geo::PlaneID DecodeCTP(CTP_t CTP)
Definition: DataStructs.cxx:89
float MCP_PFP_Cnt
Definition: TCTruth.h:64
Int_t ipart
Definition: Style.C:10
double PosSep(const Point3_t &pos1, const Point3_t &pos2)
Definition: PFPUtils.cxx:1625
double Pz(const int i=0) const
Definition: MCParticle.h:236
void StudyPiZeros(const HistStuff &hist)
Definition: TCTruth.cxx:415
CTP_t EncodeCTP(unsigned int cryo, unsigned int tpc, unsigned int plane)
Definition: DataStructs.h:45
std::vector< simb::MCParticle * > MCPartList
Definition: DataStructs.h:520
double Vz(const int i=0) const
Definition: MCParticle.h:227
std::array< double, 3 > Vector3_t
Definition: DataStructs.h:36
float fDang1
Definition: TCHist.h:86
ShowerStruct CreateSS(TjStuff &tjs, CTP_t inCTP, const std::vector< int > &tjl)
Definition: TCShower.cxx:4572
void StudyElectrons(const HistStuff &hist)
Definition: TCTruth.cxx:386
std::vector< evd::details::RawDigitInfo_t >::const_iterator end(RawDigitCacheDataClass const &cache)
TPCGeo const & TPC(unsigned int const tpc=0, unsigned int const cstat=0) const
Returns the specified TPC.
float fMCSMom
Definition: TCHist.h:86
2D representation of charge deposited in the TDC/wire plane
Definition: Hit.h:49
bool InsideTPC(const TjStuff &tjs, Point3_t &pos, geo::TPCID &inTPCID)
Definition: PFPUtils.cxx:2725
float ChgToMeV(float chg)
Definition: TCShower.cxx:4460
unsigned int ChannelID_t
Type representing the ID of a readout channel.
Definition: RawTypes.h:27
TPCID_t TPC
Index of the TPC within its cryostat.
Definition: geo_types.h:203
const simb::MCParticle * TrackIdToMotherParticle_P(int const &id)
void PrintResults(int eventNum) const
Definition: TCTruth.cxx:1083
size_type size() const
Definition: ParticleList.h:313
std::vector< int > GetAssns(const TjStuff &tjs, std::string type1Name, int id, std::string type2Name)
Definition: Utils.cxx:4079
double Vy(const int i=0) const
Definition: MCParticle.h:226
std::array< unsigned short, 3 > TruVxCounts
Definition: TCTruth.h:78
unsigned int GetMCPartListIndex(const TrajPoint &tp)
Definition: TCTruth.cxx:1499
Cosmic rays.
Definition: MCTruth.h:22
void MakeTruTrajPoint(unsigned int MCParticleListIndex, TrajPoint &tp)
Definition: TCTruth.cxx:1174
const art::Ptr< simb::MCTruth > & TrackIdToMCTruth_P(int const &id)
Beam neutrinos.
Definition: MCTruth.h:21
float Prim_TSum
Definition: TCTruth.h:65