45 if(plist.
empty())
return;
48 int sourcePtclTrackID = -1;
59 sourcePtclTrackID = trackID;
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";
71 sourcePtclTrackID = trackID;
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";
83 sourcePtclTrackID = trackID;
88 if(sourcePtclTrackID == -1) {
89 if(
tjs.
MatchTruth[1] > 0) std::cout<<
"MatchTrueHits: SourcePtcl not found\n";
96 myprt<<
"Displaying all neutrino origin MCParticles with T > 50 MeV\n";
97 myprt<<
" trackID PDGCode Mother T(MeV) ________dir_______ Process\n";
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;
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();
124 std::vector<bool> select(plist.
size(),
false);
126 unsigned int indx = 0;
130 if(theTruth->
Origin() != sourceOrigin)
continue;
136 std::vector<int> gtid(
tjs.
fHits.size(), 0);
138 for(
unsigned int iht = 0; iht <
tjs.
fHits.size(); ++iht) {
139 std::vector<sim::TrackIDE> ides;
142 geo::PlaneID planeID =
geo::PlaneID(tcHit.ArtPtr->WireID().Cryostat, tcHit.ArtPtr->WireID().TPC, tcHit.ArtPtr->WireID().Plane);
144 tcHit.StartTick, tcHit.EndTick,
145 tcHit.PeakTime, tcHit.SigmaPeakTime,
147 tcHit.PeakAmplitude, tcHit.SigmaPeakAmp,
148 tcHit.Integral, tcHit.Integral, tcHit.SigmaIntegral,
149 tcHit.Multiplicity, tcHit.LocalIndex,
150 tcHit.GoodnessOfFit, tcHit.NDOF,
153 tcHit.ArtPtr->WireID());
158 if(ides.empty())
continue;
160 for(
auto ide : ides) energy += ide.energy;
161 if(energy == 0)
continue;
165 for(
auto ide : ides) {
166 if(ide.energy > energy) {
167 hitTruTrkID = ide.trackID;
171 if(hitTruTrkID == 0)
continue;
174 if(theTruth->
Origin() != sourceOrigin)
continue;
175 unsigned int indx = 0;
179 if(mcp->
TrackId() != hitTruTrkID)
continue;
180 select[indx - 1] =
true;
184 if(!momMCP)
continue;
186 unsigned int mindx = 0;
190 select[mindx] =
true;
210 std::cout<<
"MatchTrueHits: Crazy large number of MCParticles "<<
tjs.
MCPartList.size()<<
". Ignoring this event\n";
216 unsigned int nMatch = 0;
217 for(
unsigned int iht = 0; iht <
tjs.
fHits.size(); ++iht) {
218 if(gtid[iht] == 0)
continue;
220 for(
unsigned int indx = 0; indx <
tjs.
MCPartList.size(); ++indx) {
222 if(mcp->TrackId() != gtid[iht])
continue;
223 hit.MCPartListIndex = indx;
230 std::cout<<
"MatchTrueHits: MCPartList size "<<
tjs.
MCPartList.size()<<
" nMatched hits "<<nMatch<<
"\n";
242 Point3_t primVx {{-666.0, -666.0, -666.0}};
245 primVx[0] = primMCP->Vx();
246 primVx[1] = primMCP->Vy();
247 primVx[2] = primMCP->Vz();
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();
258 std::string fcnLabel =
"SSP";
265 if(abs(mcp->PdgCode()) != 11 && abs(mcp->PdgCode()) != 111)
continue;
268 if(mcp->TrackId() != eveID)
continue;
270 auto ss3 = mcpu.MakeCheatShower(
part, primVx, truPFP);
271 if(ss3.ID == 0)
continue;
272 if(truPFP == 0)
continue;
274 std::cout<<
"Failed to store 3S"<<ss3.ID<<
"\n";
280 if(pfp.TPCID != ss3.TPCID)
continue;
282 if(pfp.PDGCode == 12 || pfp.PDGCode == 14)
continue;
284 if(pfp.PDGCode == 1111)
continue;
286 float minEnergy = 1E6;
287 for(
auto tid : pfp.TjIDs) {
291 if(energy < minEnergy) minEnergy =
energy;
293 pfpEnergy -= minEnergy;
294 pfpEnergy /= (float)(pfp.TjIDs.size() - 1);
296 unsigned short pEnd =
FarEnd(
tjs, pfp, ss3.ChgPos);
299 float costh1 = std::abs(
DotProd(pToS, ss3.Dir));
300 float costh2 =
DotProd(pToS, pfp.Dir[pEnd]);
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);
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;
317 hist.
fSep = sqrt(distToChgPos2);
322 hist.
fDang1 = acos(costh1);
323 hist.
fDang2 = acos(costh2);
329 for(
unsigned short plane = 0; plane <
tjs.
NumPlanes; ++plane) {
330 CTP_t inCTP =
EncodeCTP(ss3.TPCID.Cryostat, ss3.TPCID.TPC, plane);
332 for(
auto cid : ss3.CotIDs) {
334 if(
ss.CTP != inCTP)
continue;
338 if(ssid == 0)
continue;
342 float sep =
PosSep(tpFrom.Pos, stp1.Pos);
343 float toPos = tpFrom.Pos[0] + 0.5 * tpFrom.Dir[0] * sep;
349 if(totSep > 0) hist.
fChgFrac = chgFrac / totSep;
350 hist.
fAlong = alongTrans[0];
351 hist.
fTrans = alongTrans[1];
354 if(pfp.ID == truPFP && isBad) {
356 myprt<<
"SSP: 3S"<<ss3.ID<<
" shEnergy "<<(int)ss3Energy<<
" P"<<pfp.ID<<
" pfpEnergy "<<(
int)pfpEnergy;
362 if(pfp.ID == truPFP) {
373 if(ss3.ID == 0)
continue;
374 if(!ss3.Cheat)
continue;
375 for(
auto cid : ss3.CotIDs) {
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;
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;
409 float len =
PosSep(tj.Pts[tj.EndPt[0]].Pos, tj.Pts[tj.EndPt[1]].Pos);
421 unsigned short cnt = 0;
426 if(cnt == 1 && pdg != 111) {
427 std::cout<<
"First MC particle isn't a pizero\n";
430 if(cnt == 1)
continue;
432 double photE = 1000 * p->
E();
433 if(photE < 50)
continue;
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;
444 if(eveID != p->
TrackId())
continue;
445 unsigned short plane =
tjs.
fHits[iht].ArtPtr->WireID().Plane;
446 chgSum[plane] +=
tjs.
fHits[iht].Integral;
448 for(
unsigned short plane = 0; plane < 3; ++plane) {
449 if(chgSum[plane] == 0)
continue;
450 float chgCal = photE / chgSum[plane];
466 for(
auto& pfp :
tjs.
pfps) pfp.EffPur = 0;
470 auto const* SCE = lar::providerFrom<spacecharge::SpaceChargeService>();
474 bool neutrinoVxReconstructable =
false;
475 bool vxReconstructedNearNuVx =
false;
476 bool neutrinoPFPCorrect =
false;
478 Point3_t primVx {{-666.0, -666.0, -666.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";
492 if(
tjs.
MatchTruth[1] > 0) std::cout<<
"Found a primary particle but it is not inside any TPC\n";
495 neutrinoVxReconstructable =
true;
501 std::vector<unsigned int> mcpSelect;
503 std::vector<unsigned int> primMCPs;
507 int pdg = abs(mcp->PdgCode());
508 bool isCharged = (pdg == 11) || (pdg == 13) || (pdg == 211) || (pdg == 321) || (pdg == 2212);
509 if(!isCharged)
continue;
512 mcpSelect.push_back(
part);
514 if(mcp->Mother() != 0)
continue;
516 primMCPs.push_back(
part);
518 if(mcpSelect.empty())
return;
525 std::vector<std::pair<unsigned int, unsigned int>> moda;
526 for(
unsigned int part = 0;
part < mcpSelect.size(); ++
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);
535 for(
unsigned short dpart =
part + 1; dpart < mcpSelect.size(); ++dpart) {
537 if(dmcp->TrackId() == dtrTrackId) {
546 if(ndtr != 1)
continue;
549 if(dtr->PdgCode() != mcp->PdgCode())
continue;
551 moda.push_back(std::make_pair(mcpSelect[part], mcpSelect[dtrSelIndex]));
553 mcpSelect.erase(mcpSelect.begin() + dtrSelIndex);
563 if(moda.size() > 1) std::reverse(moda.begin(), moda.end());
566 for(
auto& md : moda)
if(md.second ==
hit.MCPartListIndex)
hit.MCPartListIndex = md.first;
571 for(
unsigned int part = 0;
part < mcpSelect.size(); ++
part) {
574 if(pdgIndex > 4)
continue;
575 float TMeV = 1000 * (mcp->E() - mcp->Mass());
576 hist.
fTruT[pdgIndex]->Fill(TMeV);
582 std::vector<std::vector<unsigned short>> ntht(
tjs.
NumPlanes);
583 for(
unsigned short plane = 0; plane <
tjs.
NumPlanes; ++plane) {
589 if(tj.AlgMod[
kKilled])
continue;
592 if(planeID.
TPC != inTPCID.
TPC)
continue;
593 tj.MCPartListIndex = UINT_MAX;
594 std::vector<unsigned int> mcpIndex, cnt;
596 for(
auto iht : tjhits) {
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);
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) {
615 tmpIndex = mcpIndex[ii];
618 if(tmpIndex == UINT_MAX)
continue;
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;
628 if(tjid[plane][tmpIndex] > 0) {
629 if(maxcnt > ntht[plane][tmpIndex]) {
631 auto& mtj =
tjs.
allTraj[tjid[plane][tmpIndex] - 1];
633 mtj.MCPartListIndex = UINT_MAX;
636 tj.MCPartListIndex = tmpIndex;
637 tjid[plane][tmpIndex] = tj.ID;
638 ntht[plane][tmpIndex] = maxcnt;
645 tj.EffPur = eff * pur;
646 tj.MCPartListIndex = tmpIndex;
647 tjid[plane][tmpIndex] = tj.ID;
648 ntht[plane][tmpIndex] = maxcnt;
652 if(neutrinoVxReconstructable) {
655 unsigned short vx3ID = 0;
658 if((pfp.PDGCode == 14 || pfp.PDGCode == 12) && pfp.Vx3ID[0] > 0) {
660 auto& vx3 =
tjs.
vtx3[pfp.Vx3ID[0] - 1];
664 Point3_t vpos = {{vx3.X, vx3.Y, vx3.Z}};
665 if(
PosSep(vpos, primVx) < 1) neutrinoPFPCorrect =
true;
669 if(vx3.ID == 0)
continue;
670 if(vx3.TPCID != inTPCID)
continue;
673 Point3_t vpos = {{vx3.X, vx3.Y, vx3.Z}};
674 float sep =
PosSep(vpos, primVx);
680 if(vx3ID > 0) vxReconstructedNearNuVx =
true;
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";
692 bool doPrt = (std::find(mcpSelect.begin(), mcpSelect.end(),
ipart) != mcpSelect.end());
695 if(mcp->PdgCode() == 111) doPrt =
true;
697 int TMeV = 1000 * (mcp->E() - mcp->Mass());
698 if(mcp->PdgCode() == 22 && mcp->Process() ==
"Decay" && TMeV > 30) doPrt =
true;
701 myprt<<std::setw(7)<<
ipart;
702 myprt<<std::setw(8)<<mcp->TrackId();
703 myprt<<std::setw(6)<<mcp->PdgCode();
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);
721 for(
unsigned short xyz = 0; xyz < 3; ++xyz) myprt<<std::setw(6)<<std::setprecision(2)<<
dir[xyz];
723 for(
unsigned short plane = 0; plane <
tjs.
NumPlanes; ++plane) {
726 myprt<<std::setw(6)<<tp.Delta;
727 startWire[plane] = tp.Pos[0];
729 myprt<<std::setw(22)<<mcp->Process();
733 for(
unsigned short plane = 0; plane < TPC.
Nplanes(); ++plane) {
736 if(mcpHits.empty())
continue;
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])) {
745 myprt<<
"_"<<mcpHits.size();
777 if(tj.AlgMod[
kKilled])
continue;
778 if(tj.MCPartListIndex == UINT_MAX)
continue;
781 if(truIndex == SHRT_MAX)
continue;
784 for(
unsigned short ipt = tj.EndPt[0]; ipt <= tj.EndPt[1]; ++ipt) {
785 auto& tp = tj.Pts[ipt];
789 if(cnt > 0) frac /= cnt;
790 float TMeV = 1000 * (mcp->E() - mcp->Mass());
791 hist.
fNearTj[truIndex]->Fill(TMeV, frac);
819 constexpr
double twopi = 2 * M_PI;
820 std::array<int, 5> recoCodeList = {{0, 11, 13, 211, 2212}};
822 if(pfp.ID == 0)
continue;
824 if(pfp.PDGCode == 1111)
continue;
826 if(pfp.MCPartListIndex == UINT_MAX)
continue;
829 if(truIndex == SHRT_MAX)
continue;
831 for(recIndex = 0; recIndex < 5; ++recIndex)
if(pfp.PDGCode == recoCodeList[recIndex])
break;
832 if(recIndex > 4)
continue;
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];
846 if(
PosSep(pfp.XYZ[1], truStart) <
PosSep(pfp.XYZ[0], truStart)) {
848 for(
unsigned short xyz = 0; xyz < 3; ++xyz) recDir[xyz] *= -1;
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()}};
863 float dang =
DeltaAngle(truDir, pfp.Dir[startEnd]);
864 while(dang > M_PI) dang -= twopi;
865 while(dang < -M_PI) dang += twopi;
877 if(mcpSelect.empty())
return;
879 unsigned int tpc = inTPCID.
TPC;
880 unsigned int cstat = inTPCID.
Cryostat;
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);
895 std::vector<std::vector<unsigned int>> tjHits(
tjs.
allTraj.size());
896 for(
unsigned short itj = 0; itj <
tjs.
allTraj.size(); ++itj) {
898 if(tj.AlgMod[
kKilled])
continue;
901 tj.MCPartListIndex = UINT_MAX;
905 for(
unsigned short isel = 0; isel < mcpSelect.size(); ++isel) {
906 if(mcpHits[isel].empty())
continue;
907 unsigned int mcpIndex = mcpSelect[isel];
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) {
923 std::vector<unsigned int> mcpPlnHits;
924 for(
auto iht : mcpHits[isel]) {
926 if(
hit.ArtPtr->WireID().Plane == plane) mcpPlnHits.push_back(iht);
929 if(mcpPlnHits.size() < 2)
continue;
931 TSums[pdgIndex] += TMeV;
934 unsigned short mtj = USHRT_MAX;
936 for(
unsigned short itj = 0; itj <
tjs.
allTraj.size(); ++itj) {
938 if(tjHits[itj].empty())
continue;
941 if(tj.CTP != inCTP)
continue;
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;
950 hist.
fEff->Fill(eff);
951 hist.
fPur->Fill(pur);
954 if(tj.MCPartListIndex != UINT_MAX && ep > tj.EffPur) {
956 tj.MCPartListIndex = mcpIndex;
963 if(mtj == USHRT_MAX) {
966 hist.
fEP_T[pdgIndex]->Fill(TMeV, 0);
970 myprt<<particleName<<
" BadEP TMeV "<<(int)TMeV<<
" No matched trajectory to isel "<<isel;
971 myprt<<
" nTrue hits "<<mcpPlnHits.size();
979 if(maxEP < tj.EffPur)
continue;
981 tj.MCPartListIndex = mcpIndex;
982 EPTSums[pdgIndex] += TMeV * tj.EffPur;
983 hist.
fEP_T[pdgIndex]->Fill(TMeV, tj.EffPur);
988 myprt<<particleName<<
" BadEP: "<<std::fixed<<std::setprecision(2)<<tj.EffPur;
989 myprt<<
" mcpIndex "<<mcpIndex;
990 myprt<<
" TMeV "<<(int)TMeV<<
" MCP hits "<<mcpPlnHits.size();
1002 std::vector<std::vector<unsigned int>> pfpHits;
1004 pfpHits.resize(
tjs.
pfps.size());
1006 for(
unsigned short ipfp = 0; ipfp <
tjs.
pfps.size(); ++ipfp) {
1008 if(pfp.ID == 0)
continue;
1010 if(pfp.TPCID != inTPCID)
continue;
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());
1022 for(
unsigned short isel = 0; isel < mcpSelect.size(); ++isel) {
1023 unsigned short mpfp = USHRT_MAX;
1025 unsigned int mcpIndex = mcpSelect[isel];
1026 if(mcpHits[isel].empty())
continue;
1028 float TMeV = 1000 * (mcp->E() - mcp->Mass());
1032 bool longMCP = (pdgIndex > 0 && pdgIndex < 5 && (float)mcpHits[isel].size() >= 2 *
tjs.
MatchTruth[3]);
1034 for(
unsigned short ipfp = 0; ipfp <
tjs.
pfps.size(); ++ipfp) {
1036 if(pfp.ID == 0)
continue;
1038 if(pfp.TPCID != inTPCID)
continue;
1040 if(pfpHits[ipfp].empty())
continue;
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;
1047 if(pfp.MCPartListIndex != UINT_MAX && ep > pfp.EffPur) {
1049 pfp.MCPartListIndex = mcpIndex;
1056 if(mpfp == USHRT_MAX) {
1059 myprt<<
"BadPFP: MCParticle "<<mcpSelect[isel]<<
" w PDGCode "<<mcp->PdgCode()<<
" T "<<(int)TMeV<<
" not reconstructed.";
1060 myprt<<
" matched Tjs:";
1062 if(tj.AlgMod[
kKilled])
continue;
1063 if(tj.MCPartListIndex == mcpIndex) myprt<<
" "<<tj.ID<<
" EP "<<std::fixed<<std::setprecision(2)<<tj.EffPur;
1070 if(maxEP < pfp.EffPur)
continue;
1072 pfp.MCPartListIndex = mcpIndex;
1089 myprt<<
"Evt "<<eventNum;
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;
1103 sum +=
TSums[pdgIndex];
1107 if(sum > 0) myprt<<
" MuPiKP "<<std::fixed<<std::setprecision(2)<<sumt / sum;
1108 myprt<<
" BadEP "<<
nBadEP;
1110 float longGood = 1 - (float)nBadEP / (
float)
nLongInPln;
1111 myprt<<
" longGood "<<std::fixed<<std::setprecision(2)<<longGood;
1116 myprt<<
" MCP cnt "<<(int)
MCP_Cnt<<
" PFP "<<std::fixed<<std::setprecision(2)<<ep;
1120 myprt<<
" PrimPFP "<<std::fixed<<std::setprecision(2)<<ep;
1124 myprt<<
" longGood "<<std::fixed<<std::setprecision(2)<<longGood;
1129 myprt<<
" NuVx correct "<<std::fixed<<std::setprecision(2)<<frac;
1139 if(nDimensions < 2 || nDimensions > 3)
return false;
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];
1147 unsigned short nPlnOK = 0;
1149 for(
unsigned short plane = 0; plane <
tjs.
NumPlanes; ++plane)
if(cntInPln[plane] > 1) ++nPlnOK;
1150 if(nPlnOK < nDimensions - 1)
return false;
1159 std::vector<unsigned int> hitVec;
1162 for(
unsigned int iht = 0; iht <
tjs.
fHits.size(); ++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);
1199 int goodID = ss3.ID;
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";
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;
1221 for(
unsigned int iht = 0; iht <
tjs.
fHits.size(); ++iht) {
1224 if(
hit.ArtPtr->WireID().TPC != tpc)
continue;
1225 if(
hit.ArtPtr->WireID().Cryostat != cstat)
continue;
1226 if(
hit.Integral <= 0)
continue;
1229 if(
hit.InTraj <= 0)
continue;
1232 if(abs(shmcp->PdgCode()) != 11)
continue;
1235 if(eid != eveID)
continue;
1236 if(shmcp->TrackId() == eid) {
1238 parentHits[
hit.ArtPtr->WireID().Plane].push_back(iht);
1242 showerHits[
hit.ArtPtr->WireID().Plane].push_back(iht);
1247 if(npar > nsh)
return ss3;
1255 std::vector<int> dummy_tjlist;
1256 for(
unsigned short plane = 0; plane <
tjs.
NumPlanes; ++plane) {
1260 ss.ShPts.resize(showerHits[plane].size());
1261 for(
unsigned short iht = 0; iht < showerHits[plane].size(); ++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;
1271 std::cout<<
"Failed to update 2S"<<
ss.ID<<
"\n";
1275 std::vector<ShowerPoint> spts;
1276 for(
auto& spt :
ss.ShPts) {
1279 float tau = along / shMaxAlong;
1282 if(tau > -1 && tau < 2) {
1283 spts.push_back(spt);
1288 std::cout<<
" 2S"<<
ss.ID<<
" shpts "<<
ss.ShPts.size()<<
" new "<<spts.size()<<
"\n";
1290 ss.NeedsUpdate =
true;
1292 std::cout<<
"Failed to update 2S"<<
ss.ID<<
"\n";
1296 std::cout<<
"Failed to store 2S"<<
ss.ID<<
"\n";
1299 ss3.CotIDs.push_back(
ss.ID);
1303 std::cout<<
"SS3 Failed...\n";
1310 std::vector<std::pair<int, int>> pcnt;
1311 for(
unsigned short plane = 0; plane <
tjs.
NumPlanes; ++plane) {
1312 for(
auto iht : parentHits[plane]) {
1314 if(
hit.InTraj <= 0)
continue;
1316 if(!tj.AlgMod[
kMat3D])
continue;
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));
1325 ++pcnt[indx].second;
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;
1337 truParentPFP = pfp.ID;
1340 std::cout<<
"truParent P"<<truParentPFP<<
" close "<<std::fixed<<std::setprecision(2)<<close<<
"\n";
1354 if(abs(mcp->PdgCode()) != 11)
continue;
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()}};
1366 energy = 1000 * (mcp->E() - mcp->Mass());
1384 if(abs(mcp->PdgCode()) != 11)
continue;
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();
1395 std::vector<int> pfplist;
1396 for(
unsigned short plane = 0; plane <
tjs.
NumPlanes; ++plane) {
1398 std::vector<unsigned int> mcphits = tm.PutMCPHitsInVector(
part, inCTP);
1399 if(mcphits.empty())
continue;
1400 for(
auto iht : mcphits) {
1402 if(
hit.InTraj <= 0)
continue;
1405 if(!tj.AlgMod[
kMat3D])
continue;
1409 if(shared.size() < 10)
continue;
1412 if(TInP.size() != 1)
continue;
1414 if(std::find(pfplist.begin(), pfplist.end(), pid) == pfplist.end()) pfplist.push_back(pid);
1417 if(pfplist.empty())
return 0;
1423 for(
auto pid : pfplist) {
1425 unsigned short nearEnd = 1 -
FarEnd(
tjs, pfp, start);
1426 float sep =
PosSep2(pfp.XYZ[nearEnd], start);
1427 if(sep > close)
continue;
1445 if(abs(mcp->PdgCode()) != 11)
continue;
1447 if(mcp->TrackId() != eveID)
continue;
1448 return MCParticleStartTjID(
part, inCTP);
1464 std::vector<std::array<int, 2>> t_cnt;
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()) {
1475 std::array<int, 2>
tmp;
1476 tmp[0] =
hit.InTraj;
1478 t_cnt.push_back(tmp);
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;
1494 return t_cnt[occIndx][0];
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()) {
1514 mcpIndex.push_back(mcpi);
1515 mcpCnt.push_back(1);
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];
1540 if(ss.
TjIDs.empty())
return UINT_MAX;
1544 for(
auto& tjid : ss.
TjIDs) {
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];
1552 ++pListCnt[
tjs.
fHits[iht].MCPartListIndex];
1557 unsigned int pIndex = UINT_MAX;
1559 for(
unsigned short ii = 0; ii < pListCnt.size(); ++ii) {
1560 if(pListCnt[ii] > nTruHits) {
1561 nTruHits = pListCnt[ii];
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];
1587 ++pListCnt[
tjs.
fHits[iht].MCPartListIndex];
1591 unsigned int pIndex = UINT_MAX;
1593 for(
unsigned short ii = 0; ii < pListCnt.size(); ++ii) {
1594 if(pListCnt[ii] > nTruHits) {
1595 nTruHits = pListCnt[ii];
double E(const int i=0) const
ROOT::Math::DisplacementVector3D< ROOT::Math::Cartesian3D< double >, ROOT::Math::GlobalCoordinateSystemTag > Vector_t
Type for representation of momenta in 3D space.
MaybeLogger_< ELseverityLevel::ELsev_info, true > LogVerbatim
bool SelectEvent
number of points to find AveChg
unsigned short nLongInPln
double Py(const int i=0) const
const std::vector< std::string > AlgBitNames
CTP_t CTP
Cryostat, TPC, Plane code.
double InShowerProbLong(double showerEnergy, double along)
std::array< double, 3 > Point3_t
const std::vector< sim::TrackIDE > HitToTrackIDEs(recob::Hit const &hit)
bool PrimaryElectronStart(Point3_t &start, Vector3_t &dir, float &energy)
unsigned int Nplanes() const
Number of planes in this tpc.
std::vector< unsigned int > PutTrajHitsInVector(Trajectory const &tj, HitStatus_t hitRequest)
float ChgFracBetween(TjStuff &tjs, Point3_t pos1, Point3_t pos2, geo::TPCID tpcid)
unsigned short PDGCodeIndex(TjStuff &tjs, int PDGCode)
void FindAlongTrans(Point3_t pos1, Vector3_t dir1, Point3_t pos2, Point2_t &alongTrans)
list_type::const_iterator const_iterator
simb::Origin_t Origin() const
The data type to uniquely identify a Plane.
Geometry information for a single TPC.
std::vector< PFPStruct > pfps
enum simb::_ev_origin Origin_t
event origin types
std::vector< unsigned int > Hits
double Px(const int i=0) const
CryostatID_t Cryostat
Index of cryostat.
int PrimaryElectronPFPID(const geo::TPCID &tpcid)
std::vector< ShowerStruct3D > showers
ShowerStruct3D MakeCheatShower(unsigned int mcpIndex, Point3_t primVx, int &truParentPFP)
double DeltaAngle(const Vector3_t v1, const Vector3_t v2)
double ShowerEnergy(const ShowerStruct3D &ss3)
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
std::string Process() const
int PrimaryElectronTjID(CTP_t inCTP)
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)
std::vector< T > SetIntersection(const std::vector< T > &set1, const std::vector< T > &set2)
int MCParticleStartTjID(unsigned int MCParticleListIndex, CTP_t inCTP)
unsigned short FarEnd(const TjStuff &tjs, const PFPStruct &pfp, const Point3_t &pos)
bool CanReconstruct(unsigned int mcpIndex, unsigned short nDimensions, const geo::TPCID &tpcid)
TH2F * fElectronLike_Len[5]
TH1F * fPFPStartAngDiff[5]
std::vector< ShowerStruct > cots
unsigned int EventsProcessed
bool UpdateShower(std::string inFcnLabel, TjStuff &tjs, ShowerStruct &ss, bool prt)
std::array< float, 2 > Point2_t
const geo::GeometryCore * geom
Vector3_t PointDirection(const Point3_t p1, const Point3_t p2)
single particles thrown at the detector
std::vector< TrajPoint > Pts
Trajectory points.
std::array< float, 5 > TSums
TrajPoint MakeBareTP(TjStuff &tjs, Point3_t &pos, Vector3_t &dir, CTP_t inCTP)
std::array< float, 5 > EPTSums
std::vector< evd::details::RawDigitInfo_t >::const_iterator begin(RawDigitCacheDataClass const &cache)
std::array< short, 5 > EPCnts
std::vector< Trajectory > allTraj
vector of all trajectories in each plane
std::vector< unsigned int > PutMCPHitsInVector(unsigned int mcpIndex, CTP_t inCTP)
double PosSep2(const Point3_t &pos1, const Point3_t &pos2)
std::string PrintHit(const TCHit &hit)
bool StoreShower(std::string inFcnLabel, TjStuff &tjs, ShowerStruct3D &ss3)
float UnitsPerTick
scale factor from Tick to WSE equivalent units
The data type to uniquely identify a TPC.
PlaneID_t Plane
Index of the plane within its TPC.
void MatchTruth(const HistStuff &hist, bool fStudyMode)
unsigned short nGoodLongMCP
bool SetMag(Vector3_t &v1, double mag)
View_t View(geo::PlaneID const &pid) const
Returns the view (wire orientation) on the channels of specified TPC plane.
void MatchAndSum(const HistStuff &hist, const std::vector< unsigned int > &mcpSelect, const geo::TPCID &inTPCID)
std::vector< float > MatchTruth
Match to MC truth.
Detector simulation of raw signals on wires.
short MCSMom(const TjStuff &tjs, const std::vector< int > &tjIDs)
std::vector< Vtx3Store > vtx3
3D vertices
double Vx(const int i=0) const
ShowerStruct3D CreateSS3(TjStuff &tjs, const geo::TPCID &tpcid)
raw::ChannelID_t PlaneWireToChannel(WireID const &wireid) const
Returns the ID of the TPC channel connected to the specified wire.
double DotProd(const Vector3_t &v1, const Vector3_t &v2)
std::string to_string(Flag_t< Storage > const flag)
Convert a flag into a stream (shows its index).
std::vector< TCHit > fHits
std::string PrintHitShort(const TCHit &hit)
TProfile * fChgToMeV_Etru
geo::PlaneID DecodeCTP(CTP_t CTP)
double PosSep(const Point3_t &pos1, const Point3_t &pos2)
double Pz(const int i=0) const
void StudyPiZeros(const HistStuff &hist)
CTP_t EncodeCTP(unsigned int cryo, unsigned int tpc, unsigned int plane)
std::vector< simb::MCParticle * > MCPartList
double Vz(const int i=0) const
std::array< double, 3 > Vector3_t
ShowerStruct CreateSS(TjStuff &tjs, CTP_t inCTP, const std::vector< int > &tjl)
const sim::ParticleList & ParticleList()
void StudyElectrons(const HistStuff &hist)
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.
2D representation of charge deposited in the TDC/wire plane
bool InsideTPC(const TjStuff &tjs, Point3_t &pos, geo::TPCID &inTPCID)
float ChgToMeV(float chg)
unsigned int ChannelID_t
Type representing the ID of a readout channel.
TPCID_t TPC
Index of the TPC within its cryostat.
const simb::MCParticle * TrackIdToMotherParticle_P(int const &id)
void PrintResults(int eventNum) const
std::vector< int > GetAssns(const TjStuff &tjs, std::string type1Name, int id, std::string type2Name)
double Vy(const int i=0) const
std::array< unsigned short, 3 > TruVxCounts
unsigned int GetMCPartListIndex(const TrajPoint &tp)
void MakeTruTrajPoint(unsigned int MCParticleListIndex, TrajPoint &tp)
const art::Ptr< simb::MCTruth > & TrackIdToMCTruth_P(int const &id)