36 : fCaloAlg(pset.
get<
fhicl::ParameterSet>(
"CaloAlg")), fMVAReader(
"Silent")
40 bool badinput =
false;
47 if (pset.
has_key(
"Mode")) userMode = pset.
get<
short>(
"Mode");
51 if (pset.
has_key(
"StudyMode")) {
52 std::cout <<
"StudyMode is not valid anymore. Try Study: 1 or Study: 2, etc/n";
55 userMode = pset.
get<
short>(
"Study");
61 if (pset.
has_key(
"SaveShowerTree"))
65 std::vector<std::string> skipAlgsVec;
66 if (pset.
has_key(
"SkipAlgs")) skipAlgsVec = pset.
get<std::vector<std::string>>(
"SkipAlgs");
67 std::vector<std::string> debugConfigVec;
68 if (pset.
has_key(
"DebugConfig"))
69 debugConfigVec = pset.
get<std::vector<std::string>>(
"DebugConfig");
73 std::vector<float> aveHitRMS;
74 if (pset.
has_key(
"AveHitRMS")) aveHitRMS = pset.
get<std::vector<float>>(
"AveHitRMS");
76 if (!aveHitRMS.empty()) {
83 tcc.
minPts = pset.
get<std::vector<unsigned short>>(
"MinPts");
87 tcc.
chargeCuts = pset.
get<std::vector<float>>(
"ChargeCuts", {3, 0.15, 0.25});
97 tcc.
muonTag = pset.
get<std::vector<short>>(
"MuonTag", {-1, -1, -1, -1});
99 tcc.
showerTag = pset.
get<std::vector<float>>(
"ShowerTag", {-1, -1, -1, -1, -1, -1});
100 std::string fMVAShowerParentWeights =
"NA";
101 pset.
get_if_present<std::string>(
"MVAShowerParentWeights", fMVAShowerParentWeights);
103 if (pset.
has_key(
"MatchTruth")) {
104 std::cout <<
"MatchTruth is not used. Use ClusterAnaV2 or DebugConfig to configure\n";
106 tcc.
vtx2DCuts = pset.
get<std::vector<float>>(
"Vertex2DCuts", {-1, -1, -1, -1, -1, -1, -1});
107 tcc.
vtx3DCuts = pset.
get<std::vector<float>>(
"Vertex3DCuts", {-1, -1});
109 tcc.
match3DCuts = pset.
get<std::vector<float>>(
"Match3DCuts", {-1, -1, -1, -1, -1});
122 <<
"Bad input from fcl file. Vector lengths for MinPtsFit, MaxAngleRange and MinMCSMom " 123 "should be defined for each reconstruction pass";
127 <<
"Vertex2DCuts must be size 11\n 0 = Max length definition for short TJs\n 1 = Max " 128 "vtx-TJ sep short TJs\n 2 = Max vtx-TJ sep long TJs\n 3 = Max position pull for >2 " 129 "TJs\n 4 = Max vtx position error\n 5 = Min MCSMom for one of two TJs\n 6 = Min " 130 "fraction of wires hit btw vtx and Tjs\n 7 = Min Score\n 8 = min ChgFrac at a vtx or " 131 "merge point\n 9 = max MCSMom asymmetry, 10 = require chg on wires btw vtx and tjs in " 140 <<
"Vertex3DCuts must be size > 2\n 0 = 2D Vtx max dX (cm)\n 1 = 2D Vtx max pull\n 2 = max " 141 "3D separation (cm) btw PFP and vertex";
143 std::cout <<
"WARNING: Vertex3DCuts is size 2 but should be size 3, where Vertex3DCuts[2] = " 144 "max 3D separation (cm) btw a PFP and a 3D vertex. Setting it to 3 cm\n";
149 <<
"KinkCuts must be size 3\n 0 = Number of points to fit at the end of the trajectory\n 1 " 150 "= Minimum kink significance\n 2 = Use charge in significance calculation? (yes if > " 151 "0)\n 3 = 3D kink fit length (cm). \nYou are using an out-of-date specification?\n";
157 <<
"Are you using an out-of-date specification for KinkCuts? KinkCuts[0] is the number of " 163 <<
"ChargeCuts must be size 3\n 0 = Charge pull cut\n 1 = Min allowed fractional chg RMS\n " 164 "2 = Max allowed fractional chg RMS";
168 <<
"MuonTag must be size 4\n 0 = minPtsFit\n 1 = minMCSMom\n 2= maxWireSkipNoSignal\n 3 = " 169 "min delta ray length for tagging\n 4 = dress muon window size (optional)";
172 <<
"DeltaRayTag must be size 3\n 0 = Max endpoint sep\n 1 = min MCSMom\n 2 = max MCSMom";
175 <<
"ChkStopCuts must be size 3\n 0 = Min Charge ratio\n 1 = Charge slope pull cut\n 2 = " 176 "Charge fit chisq cut\n 3 = BraggSplit FOM (optional)";
179 <<
"ShowerTag must be size 13\n 0 = Mode\n 1 = max MCSMom\n 2 = max separation (WSE " 180 "units)\n 3 = Max angle diff\n 4 = Factor * rms width\n 5 = Min half width\n 6 = min " 181 "total Tps\n 7 = Min Tjs\n 8 = max parent FOM\n 9 = max direction FOM 10 = max " 182 "AspectRatio\n 11 = min Score to preserve a vertex\n 12 = Debug showers in CTP\n";
183 std::cout <<
" Fixing this problem...";
192 <<
"Match3DCuts must be size 5\n 0 = dx (cm) matching cut\n 1 = max number of 3D " 193 "combinations\n 2 = min length for 2-view match\n 3 = number of TP3Ds in each plane to " 194 "fit in each PFP section\n 4 = max pull for accepting TP3Ds in sections\n 5 = max " 195 "ChiDOF for a SectionFit";
199 mf::LogVerbatim(
"TC") <<
"Last element of AngleRange != 90 degrees. Fixing it\n";
215 for (
auto strng : debugConfigVec) {
219 for (
unsigned short ib = 0; ib <
AlgBitNames.size(); ++ib) {
228 std::cout <<
"DecodeDebugString failed: " << strng <<
"\n";
236 if (range < 0 || range > 90)
238 <<
"Invalid angle range " << range <<
" Must be 0 - 90 degrees";
248 <<
"Increase the size of UseAlgs to at least " <<
kAlgBitSize;
257 <<
"Increase the size of EndFlag to at least " <<
kFlagBitSize;
259 bool printHelp =
false;
260 for (
unsigned short ib = 0; ib <
AlgBitNames.size(); ++ib)
267 for (
auto strng : skipAlgsVec) {
269 if (strng ==
"All") {
271 for (
unsigned short ib = 0; ib <
AlgBitNames.size(); ++ib)
276 for (
unsigned short ib = 0; ib <
AlgBitNames.size(); ++ib) {
284 std::cout <<
"******* Unknown SkipAlgs input string '" << strng <<
"'\n";
289 std::cout <<
"Valid AlgNames:";
291 std::cout <<
" " << strng;
293 std::cout <<
"Or specify All to turn all algs off\n";
296 if (fMVAShowerParentWeights !=
"NA" &&
tcc.
showerTag[0] > 0)
316 tcc.
geom = lar::providerFrom<geo::Geometry>();
338 thr = {UINT_MAX, UINT_MAX};
341 unsigned int tpc =
hit.WireID().TPC;
351 std::vector<unsigned int>& hitsInSlice,
357 if (hitsInSlice.size() < 2)
return;
360 if (!
CreateSlice(clockData, detProp, hitsInSlice, sliceID))
return;
373 <<
" AveHitRMS vector size != the number of planes ";
375 std::cout <<
"Reconstruct " << hitsInSlice.size() <<
" hits in Slice " << sliceID
376 <<
" in TPC " << slc.TPCID.TPC <<
"\n";
377 for (
unsigned short plane = 0; plane < slc.nPlanes; ++plane) {
378 CTP_t inCTP =
EncodeCTP(slc.TPCID.Cryostat, slc.TPCID.TPC, plane);
380 if (!slc.isValid)
return;
389 for (
unsigned short plane = 0; plane < slc.nPlanes; ++plane) {
390 CTP_t inCTP =
EncodeCTP(slc.TPCID.Cryostat, slc.TPCID.TPC, plane);
392 std::cout <<
"RTC: ChkVtxAssociations found an error\n";
405 std::cout <<
"SHOWER TREE STAGE NUM SIZE: " <<
stv.
StageNum.size() << std::endl;
411 mf::LogVerbatim(
"TC") <<
"RunTrajCluster failed in MakeAllTrajClusters";
421 for (
auto& tj : slc.tjs) {
422 for (
unsigned short ib = 0; ib <
AlgBitNames.size(); ++ib)
427 slc.mallTraj.resize(0);
441 if (nwires > slc.
nWires[plane])
return;
447 for (
unsigned short pass = 0; pass <
tcc.
minPtsFit.size(); ++pass) {
448 for (
unsigned int ii = 0; ii < nwires; ++ii) {
450 unsigned int iwire = 0;
451 unsigned int jwire = 0;
459 iwire = slc.
lastWire[plane] - ii - 1;
462 if (iwire > slc.
wireHitRange[plane].size() - 1)
continue;
463 if (jwire > slc.
wireHitRange[plane].size() - 1)
continue;
465 if (slc.
wireHitRange[plane][iwire].first == UINT_MAX)
continue;
466 if (slc.
wireHitRange[plane][jwire].first == UINT_MAX)
continue;
467 unsigned int ifirsthit = slc.
wireHitRange[plane][iwire].first;
468 unsigned int ilasthit = slc.
wireHitRange[plane][iwire].second;
469 unsigned int jfirsthit = slc.
wireHitRange[plane][jwire].first;
470 unsigned int jlasthit = slc.
wireHitRange[plane][jwire].second;
471 if (ifirsthit > slc.
slHits.size() || ilasthit > slc.
slHits.size()) {
472 std::cout <<
"RAT: bad hit range " << ifirsthit <<
" " << ilasthit <<
" size " 473 << slc.
slHits.size() <<
" inCTP " << inCTP <<
"\n";
476 for (
unsigned int iht = ifirsthit; iht <= ilasthit; ++iht) {
479 mf::LogVerbatim(
"TC") <<
"+++++++ Pass " << pass <<
" Found debug hit " 484 if (slc.
slHits[iht].InTraj != 0)
continue;
488 std::cout <<
"RAT: Bad allHitsIndex\n";
493 unsigned int fromWire = iHit.WireID().Wire;
494 float fromTick = iHit.PeakTime();
495 float iqtot = iHit.Integral();
496 float hitsRMSTick = iHit.RMS();
497 std::vector<unsigned int> iHitsInMultiplet;
500 iHitsInMultiplet.resize(1);
501 iHitsInMultiplet[0] = iht;
505 if (iHitsInMultiplet.empty())
continue;
507 if (iHitsInMultiplet.size() > 4)
continue;
509 if (iHitsInMultiplet.size() > 1) {
513 if (hitsRMSTick == 0)
continue;
514 bool fatIHit = (hitsRMSTick > maxHitsRMS);
516 mf::LogVerbatim(
"TC") <<
" hit RMS " << iHit.RMS() <<
" BB Multiplicity " 517 << iHitsInMultiplet.size() <<
" AveHitRMS[" << plane <<
"] " 518 <<
evt.
aveHitRMS[plane] <<
" HitsRMSTick " << hitsRMSTick
519 <<
" fatIHit " << fatIHit;
520 for (
unsigned int jht = jfirsthit; jht <= jlasthit; ++jht) {
522 if (slc.
slHits[iht].InTraj != 0)
break;
523 if (slc.
slHits[jht].InTraj != 0)
continue;
525 for (
auto& slHit : slc.
slHits) {
526 if (slHit.InTraj < 0) {
527 std::cout <<
"RAT: Dirty hit " <<
PrintHit(slHit) <<
" EventsProcessed " 532 unsigned int toWire = jwire;
535 float toTick = jHit.PeakTime();
536 float jqtot = jHit.Integral();
537 std::vector<unsigned int> jHitsInMultiplet;
540 jHitsInMultiplet.resize(1);
541 jHitsInMultiplet[0] = jht;
545 if (jHitsInMultiplet.empty())
continue;
547 if (jHitsInMultiplet.size() > 4)
continue;
550 if (hitsRMSTick == 0)
continue;
551 bool fatJHit = (hitsRMSTick > maxHitsRMS);
554 if ((fatIHit && !fatJHit) || (!fatIHit && fatJHit)) {
continue; }
558 if (jHitsInMultiplet.size() > 1)
561 bool hitsOK =
TrajHitsOK(slc, iHitsInMultiplet, jHitsInMultiplet);
563 if (!hitsOK)
continue;
566 if (!
StartTraj(work, fromWire, fromTick, toWire, toTick, inCTP, pass))
continue;
569 std::cout <<
"RAT: StartTraj major failure\n";
572 if (work.
Pts.empty()) {
585 std::cout <<
"RAT: AddHits major failure\n";
588 if (!sigOK || work.
Pts[0].Chg == 0) {
594 work.
Pts[0].Pos = work.
Pts[0].HitPos;
603 std::cout <<
"RAT: StepAway major failure\n";
612 std::cout <<
"RAT: CheckTraj major failure\n";
617 <<
"-" << work.
EndPt[1] <<
" IsGood " << work.
IsGood;
620 <<
"StepAway done: IsGood " << work.
IsGood <<
" NumPtsWithCharge " 633 auto& tj = slc.
tjs[slc.
tjs.size() - 1];
636 std::cout <<
"RAT: InTrajOK major failure. " << tj.ID <<
"\n";
649 for (
auto tp :
seeds) {
650 unsigned short nAvailable = 0;
651 for (
unsigned short ii = 0; ii < tp.Hits.size(); ++ii) {
652 if (!tp.UseHit[ii])
continue;
653 unsigned int iht = tp.Hits[ii];
654 if (slc.
slHits[iht].InTraj == 0) ++nAvailable;
661 if (nAvailable == 0)
continue;
664 for (
unsigned short ii = 0; ii < tp.Hits.size(); ++ii) {
665 if (!tp.UseHit[ii])
continue;
666 unsigned int iht = tp.Hits[ii];
667 if (slc.
slHits[iht].InTraj == 0) slc.
slHits[iht].InTraj = work.
ID;
671 if (tp.Dir[0] < 0) work.
StepDir = -1;
678 work.
Pts.push_back(tp);
683 std::cout <<
"RAT: CheckTraj major failure\n";
697 auto& tj = slc.
tjs[slc.
tjs.size() - 1];
722 for (
auto& tj : slc.
tjs) {
723 if (tj.AlgMod[
kKilled])
continue;
724 if (tj.PDGCode != 13)
continue;
739 std::cout <<
"RAT: MakeJunkVertices major failure\n";
744 for (
unsigned short ivx = 0; ivx < slc.
vtxs.size(); ++ivx) {
745 auto& vx2 = slc.
vtxs[ivx];
746 if (vx2.ID == 0)
continue;
747 if (vx2.CTP != inCTP)
continue;
759 std::cout <<
"RAT: ChkVtxAssociations found an error. Events processed " 776 for (
auto& slHit : slc.
slHits)
777 if (slHit.InTraj < 0) slHit.InTraj = 0;
781 std::vector<unsigned int> tHits;
783 for (
unsigned int iwire = slc.
firstWire[plane]; iwire < slc.
lastWire[plane] - 3; ++iwire) {
787 unsigned int jwire = iwire + 1;
790 unsigned int ifirsthit = slc.
wireHitRange[plane][iwire].first;
791 unsigned int ilasthit = slc.
wireHitRange[plane][iwire].second;
792 unsigned int jfirsthit = slc.
wireHitRange[plane][jwire].first;
793 unsigned int jlasthit = slc.
wireHitRange[plane][jwire].second;
794 for (
unsigned int iht = ifirsthit; iht <= ilasthit; ++iht) {
795 if (iht >= slc.
slHits.size())
break;
796 auto& islHit = slc.
slHits[iht];
797 if (islHit.InTraj != 0)
continue;
798 std::vector<unsigned int> iHits;
802 if (prt)
mf::LogVerbatim(
"TC") <<
"FJT: debug iht multiplet size " << iHits.size();
803 if (iHits.empty())
continue;
804 for (
unsigned int jht = jfirsthit; jht <= jlasthit; ++jht) {
805 if (jht >= slc.
slHits.size())
break;
806 auto& jslHit = slc.
slHits[jht];
807 if (jslHit.InTraj != 0)
continue;
808 if (prt &&
HitSep2(slc, iht, jht) < 100)
812 std::vector<unsigned int> jHits;
814 if (jHits.empty())
continue;
819 for (
auto iht : iHits)
820 if (slc.
slHits[iht].InTraj == 0) tHits.push_back(iht);
821 for (
auto jht : jHits)
822 if (slc.
slHits[jht].InTraj == 0) tHits.push_back(jht);
823 for (
auto tht : tHits)
824 slc.
slHits[tht].InTraj = -4;
826 if (iwire != 0) { loWire = iwire - 1; }
830 unsigned int hiWire = jwire + 1;
831 if (hiWire > slc.
nWires[plane])
break;
832 unsigned short nit = 0;
834 bool hitsAdded =
false;
835 for (
unsigned int kwire = loWire; kwire <= hiWire; ++kwire) {
836 if (slc.
wireHitRange[plane][kwire].first == UINT_MAX)
continue;
837 if (slc.
wireHitRange[plane][kwire].second == UINT_MAX)
continue;
838 unsigned int kfirsthit = slc.
wireHitRange[plane][kwire].first;
839 unsigned int klasthit = slc.
wireHitRange[plane][kwire].second;
840 for (
unsigned int kht = kfirsthit; kht <= klasthit; ++kht) {
841 if (kht >= slc.
slHits.size())
continue;
842 if (slc.
slHits[kht].InTraj != 0)
continue;
844 if (std::find(tHits.begin(), tHits.end(), kht) != tHits.end())
continue;
849 for (
auto jht : jHits) {
850 if (slc.
slHits[jht].InTraj != 0)
continue;
851 tHits.push_back(jht);
852 slc.
slHits[jht].InTraj = -4;
853 if (kwire > hiWire) hiWire = kwire;
854 if (kwire < loWire) loWire = kwire;
859 if (!hitsAdded)
break;
862 if (hiWire >= slc.
nWires[plane])
break;
865 for (
auto iht : tHits)
866 slc.
slHits[iht].InTraj = 0;
867 if (tHits.size() < 3)
continue;
870 myprt <<
"FJT: tHits";
871 for (
auto tht : tHits)
876 if (
IsGhost(slc, tHits))
break;
881 if (slc.
slHits[jht].InTraj > 0)
break;
898 unsigned short itj = 0;
899 std::vector<unsigned int> tHits;
900 std::vector<unsigned int> atHits;
901 for (
auto& tj : slc.
tjs) {
903 if (tj.AlgMod[
kKilled])
continue;
905 for (
auto& tp : tj.Pts) {
906 if (tp.Hits.size() > 16) {
909 <<
"ChkInTraj: More than 16 hits created a UseHit bitset overflow\n";
911 std::cout <<
"ChkInTraj major failure\n";
916 std::cout << someText <<
" ChkInTraj hit size mis-match in tj ID " << tj.ID
918 for (
unsigned short ib = 0; ib <
AlgBitNames.size(); ++ib)
919 if (tj.AlgMod[ib]) std::cout <<
" " <<
AlgBitNames[ib];
924 if (tHits.size() < 2) {
925 mf::LogVerbatim(
"TC") << someText <<
" ChkInTraj: Insufficient hits in traj " << tj.ID;
930 std::sort(tHits.begin(), tHits.end());
932 for (iht = 0; iht < slc.
slHits.size(); ++iht) {
933 if (slc.
slHits[iht].InTraj == tID) atHits.push_back(iht);
935 if (atHits.size() < 2) {
936 mf::LogVerbatim(
"TC") << someText <<
" ChkInTraj: Insufficient hits in atHits in traj " 937 << tj.ID <<
" Killing it";
941 if (!std::equal(tHits.begin(), tHits.end(), atHits.begin())) {
943 myprt << someText <<
" ChkInTraj failed: inTraj - UseHit mis-match for tj ID " << tID
944 <<
" tj.WorkID " << tj.WorkID <<
" atHits size " << atHits.size() <<
" tHits size " 945 << tHits.size() <<
" in CTP " << tj.CTP <<
"\n";
946 myprt <<
"AlgMods: ";
947 for (
unsigned short ib = 0; ib <
AlgBitNames.size(); ++ib)
948 if (tj.AlgMod[ib]) myprt <<
" " <<
AlgBitNames[ib];
950 myprt <<
"index inTraj UseHit \n";
951 for (iht = 0; iht < atHits.size(); ++iht) {
952 myprt <<
"iht " << iht <<
" " <<
PrintHit(slc.
slHits[atHits[iht]]);
953 if (iht < tHits.size()) myprt <<
" " <<
PrintHit(slc.
slHits[tHits[iht]]);
954 if (atHits[iht] != tHits[iht]) myprt <<
" <<< " << atHits[iht] <<
" != " << tHits[iht];
958 if (tHits.size() > atHits.size()) {
959 for (iht = atHits.size(); iht < atHits.size(); ++iht) {
960 myprt <<
"atHits " << iht <<
" " <<
PrintHit(slc.
slHits[atHits[iht]]) <<
"\n";
966 for (
unsigned short end = 0;
end < 2; ++
end) {
967 if (tj.VtxID[
end] > slc.
vtxs.size()) {
968 mf::LogVerbatim(
"TC") << someText <<
" ChkInTraj: Bad VtxID " << tj.ID;
969 std::cout << someText <<
" ChkInTraj: Bad VtxID " << tj.ID <<
" vtx size " 970 << slc.
vtxs.size() <<
"\n";
985 std::vector<recob::Hit>& newHitCol,
986 std::vector<unsigned int>& newHitAssns)
const 991 if (tpHits.empty())
return;
994 if (tpHits.size() == 1) {
996 std::cout <<
"MergeTPHits Bad input hit index " << tpHits[0] <<
" allHits size " 1001 newHitAssns[tpHits[0]] = newHitCol.size() - 1;
1006 std::vector<unsigned int> wires;
1007 std::vector<std::vector<unsigned int>> wireHits;
1009 wires.push_back(firstHit.WireID().Wire);
1010 std::vector<unsigned int>
tmp(1, tpHits[0]);
1011 wireHits.push_back(tmp);
1012 for (
unsigned short ii = 1; ii < tpHits.size(); ++ii) {
1014 unsigned int wire =
hit.WireID().Wire;
1015 unsigned short indx = 0;
1016 for (indx = 0; indx < wires.size(); ++indx)
1017 if (wires[indx] == wire)
break;
1018 if (indx == wires.size()) {
1019 wires.push_back(wire);
1020 wireHits.resize(wireHits.size() + 1);
1022 wireHits[indx].push_back(tpHits[ii]);
1026 for (
unsigned short indx = 0; indx < wireHits.size(); ++indx) {
1027 auto& hitsOnWire = wireHits[indx];
1029 for (
unsigned short ii = 0; ii < hitsOnWire.size(); ++ii) {
1030 newHitAssns[hitsOnWire[ii]] = newHitCol.size() - 1;
1046 if (tpHits.size() == 1) {
1048 std::cout <<
"MergeTPHits Bad input hit index " << tpHits[0] <<
" allHits size " 1060 oldHit.SigmaPeakTime(),
1062 oldHit.PeakAmplitude(),
1063 oldHit.SigmaPeakAmplitude(),
1066 oldHit.SigmaIntegral(),
1072 oldHit.SignalType(),
1076 double integral = 0;
1077 double sIntegral = 0;
1078 double peakTime = 0;
1079 double sPeakTime = 0;
1081 double sPeakAmp = 0;
1085 for (
auto allHitsIndex : tpHits) {
1088 if (
hit.StartTick() < startTick) startTick =
hit.StartTick();
1089 if (
hit.EndTick() > endTick) endTick =
hit.EndTick();
1090 double intgrl =
hit.Integral();
1091 sPeakTime += intgrl *
hit.SigmaPeakTime();
1092 sPeakAmp += intgrl *
hit.SigmaPeakAmplitude();
1093 sumADC +=
hit.SummedADC();
1095 sIntegral += intgrl *
hit.SigmaIntegral();
1098 if (integral <= 0) {
1099 std::cout <<
"MergeTPHits found bad hit integral " << integral <<
"\n";
1104 std::vector<double> shape(endTick - startTick + 1, 0.);
1105 for (
auto allHitsIndex : tpHits) {
1107 double peakTick =
hit.PeakTime();
1108 double rms =
hit.RMS();
1109 double peakAmp =
hit.PeakAmplitude();
1111 double arg = ((double)
tick - peakTick) / rms;
1112 unsigned short indx =
tick - startTick;
1113 shape[indx] += peakAmp * exp(-0.5 * arg * arg);
1121 unsigned short indx =
tick - startTick;
1122 sigsum += shape[indx];
1123 sigsumt += shape[indx] *
tick;
1126 peakTime = sigsumt / sigsum;
1128 sPeakTime /= integral;
1133 double dTick =
tick - peakTime;
1134 unsigned short indx =
tick - startTick;
1135 sigsumt += shape[indx] * dTick * dTick;
1137 double rms = std::sqrt(sigsumt / sigsum);
1140 double chgNorm = 2.507 * firstHit.PeakAmplitude() * firstHit.RMS() / firstHit.Integral();
1142 peakAmp = (float)(integral * chgNorm / (2.507 * rms));
1144 sPeakAmp /= integral;
1146 startTick = peakTime - 3 * rms;
1147 endTick = peakTime + 3 * rms;
1166 firstHit.SignalType(),
1215 std::vector<unsigned int>& hitsInSlice,
1221 if (hitsInSlice.size() < 2)
return false;
1225 slc.
slHits.resize(hitsInSlice.size());
1227 unsigned int cstat = 0;
1228 unsigned int tpc = UINT_MAX;
1229 unsigned int cnt = 0;
1230 std::vector<unsigned int> nHitsInPln;
1231 for (
auto iht : hitsInSlice) {
1232 if (iht >= (*
evt.
allHits).size())
return false;
1235 cstat =
hit.WireID().Cryostat;
1236 tpc =
hit.WireID().TPC;
1241 if (
hit.WireID().Cryostat != cstat ||
hit.WireID().TPC != tpc)
return false;
1242 slc.
slHits[cnt].allHitsIndex = iht;
1243 ++nHitsInPln[
hit.WireID().Plane];
1247 for (
auto hip : nHitsInPln)
1248 if (hip < 2)
return false;
1258 if (
tcc.
dbgSlc) std::cout <<
"Enabled debugging in sub-slice " <<
slices.size() - 1 <<
"\n";
1266 for (
unsigned int iht = 0; iht < slc.
slHits.size(); ++iht) {
1267 unsigned int ahi = slc.
slHits[iht].allHitsIndex;
1269 std::cout <<
"CreateSlice: Bad allHitsIndex " << ahi <<
" " << (*
evt.
allHits).
size()
1284 for (
auto& slc :
slices) {
1285 if (!slc.isValid)
continue;
1287 for (
auto& pfp : slc.pfps) {
1288 if (pfp.ID <= 0)
continue;
1289 pfp.TjUIDs.resize(pfp.TjIDs.size());
1290 for (
unsigned short ii = 0; ii < pfp.TjIDs.size(); ++ii) {
1292 if (pfp.TjIDs[ii] <= 0 || pfp.TjIDs[ii] > (
int)slc.tjs.size()) {
1293 std::cout <<
"FinishEvent found an invalid T" << pfp.TjIDs[ii] <<
" in P" << pfp.UID
1295 slc.isValid =
false;
1298 auto& tj = slc.tjs[pfp.TjIDs[ii] - 1];
1299 pfp.TjUIDs[ii] = tj.UID;
1306 for (
auto& slc : slices)
Expect tracks entering from the front face. Don't create neutrino PFParticles.
calo::CalorimetryAlg * caloAlg
void ClearResults()
Deletes all the results.
void CheckTrajBeginChg(TCSlice &slc, unsigned short itj)
float AveChg
Calculated using ALL hits.
std::vector< int > EnvStage
MaybeLogger_< ELseverityLevel::ELsev_info, true > LogVerbatim
std::vector< int > IsShowerParent
std::vector< Trajectory > tjs
vector of all trajectories in each plane
void StepAway(TCSlice &slc, Trajectory &tj)
void ConfigureMVA(TCConfig &tcc, std::string fMVAShowerParentWeights)
unsigned int NTPC(CryostatID const &cryoid=cryostat_zero) const
Returns the total number of TPCs in the specified cryostat.
bool TrajHitsOK(TCSlice const &slc, const std::vector< unsigned int > &iHitsInMultiplet, const std::vector< unsigned int > &jHitsInMultiplet)
std::vector< float > kinkCuts
kink finder algorithm
bool IsGhost(TCSlice &slc, Trajectory &tj)
Utilities related to art service access.
bool InTrajOK(TCSlice &slc, std::string someText)
std::vector< float > EndWir
std::vector< float > EndAng
std::vector< unsigned int > PutTrajHitsInVector(const Trajectory &tj, HitStatus_t hitRequest)
const std::vector< std::string > AlgBitNames
short recoTPC
only reconstruct in the seleted TPC
std::vector< unsigned short > minPtsFit
Reconstruct in several passes.
std::vector< float > vtx3DCuts
2D vtx -> 3D vtx matching cuts
std::vector< float > qualityCuts
Min points/wire, min consecutive pts after a gap.
void Find3DVertices(detinfo::DetectorPropertiesData const &detProp, TCSlice &slc)
std::vector< float > BeginTim
void AddHits(TCSlice &slc, Trajectory &tj, unsigned short ipt, bool &sigOK)
std::vector< float > BeginAng
void PrintTrajectory(std::string someText, const TCSlice &slc, const Trajectory &tj, unsigned short tPoint)
void SplitTrajCrossingVertices(TCSlice &slc, CTP_t inCTP)
bool StartTraj(TCSlice const &slc, Trajectory &tj, unsigned int fromhit, unsigned int tohit, unsigned short pass)
std::vector< std::vector< std::pair< unsigned int, unsigned int > > > wireHitRange
void Reconcile2Vs(TCSlice &slc)
int ParentID
ID of the parent, or the ID of the Tj this one was merged with if it is killed.
step from US -> DS (true) or DS -> US (false)
bool StoreTraj(TCSlice &slc, Trajectory &tj)
void RunTrajClusterAlg(detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, std::vector< unsigned int > &hitsInSlice, int sliceID)
calo::CalorimetryAlg fCaloAlg
float maxWireSkipWithSignal
max number of wires to skip with a signal on them
call study functions to develop cuts, etc
for(Int_t i=0;i< nentries;i++)
void FillWireHitRange(geo::TPCID inTPCID)
bool AttachAnyTrajToVertex(TCSlice &slc, unsigned short ivx, bool prt)
std::vector< std::pair< unsigned int, unsigned int > > tpcSrcHitRange
bool MakeJunkTraj(TCSlice &slc, std::vector< unsigned int > tHits)
std::vector< float > electronTag
bool FindShowers3D(detinfo::DetectorPropertiesData const &detProp, TCSlice &slc)
bool dbgSlc
debug only in the user-defined slice? default is all slices
void TagShowerLike(TCSlice &slc, const CTP_t &inCTP)
std::vector< unsigned int > lastWire
the last wire with a hit
bool CreateSlice(detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, std::vector< unsigned int > &hitsInSlice, int sliceID)
bool LongPulseHit(const recob::Hit &hit)
unsigned int eventsProcessed
bool IsGood
set false if there is a failure or the Tj fails quality cuts
std::string PrintPos(const TrajPoint &tp)
bool doForecast
See TCMode_t above.
void MakePFPTjs(TCSlice &slc)
std::vector< float > angleRanges
list of max angles for each angle range
const std::vector< std::string > EndFlagNames
std::vector< float > showerTag
shower-like trajectory tagging + shower reconstruction
unsigned short Pass
the pass on which it was created
int TDCtick_t
Type representing a TDC tick.
std::vector< float > EndTim
std::vector< int > ShowerID
decltype(auto) constexpr end(T &&obj)
ADL-aware version of std::end.
decltype(auto) constexpr size(T &&obj)
ADL-aware version of std::size.
std::vector< unsigned short > maxAngleCode
max allowed angle code for each pass
bool SetInputHits(std::vector< recob::Hit > const &inputHits, unsigned int run, unsigned int event)
void DefineShTree(TTree *t)
std::vector< unsigned int > fAlgModCount
float projectionErrFactor
unsigned int Hit
set to the hit index in evt.allHits if a Plane:Wire:Tick match is found
std::vector< std::string > StageName
float HitSep2(const TCSlice &slc, unsigned int iht, unsigned int jht)
int Cryostat
Select Cryostat.
float HitsRMSTick(const TCSlice &slc, const std::vector< unsigned int > &hitsInMultiplet, HitStatus_t hitRequest)
short nPtsAve
dump trajectory points
void PFPVertexCheck(TCSlice &slc)
bool dbgStp
debug stepping using debug.Cryostat, debug.TPC, etc
recob::Hit MergeTPHitsOnWire(std::vector< unsigned int > &tpHits) const
don't mess with this line
float unitsPerTick
scale factor from Tick to WSE equivalent units
std::vector< short > minMCSMom
Min MCSMom for each pass.
bool BraggSplit(TCSlice &slc, unsigned short itj)
std::vector< short > BeginVtx
CTP_t CTP
Cryostat, TPC, Plane code.
bool dbg2V
debug 2D vertex finding
void GetHitMultiplet(const TCSlice &slc, unsigned int theHit, std::vector< unsigned int > &hitsInMultiplet, bool useLongPulseHits)
std::vector< float > aveHitRMS
average RMS of an isolated hit
std::vector< TrajPoint > Pts
Trajectory points.
TMVA::Reader * showerParentReader
bool ChkVtxAssociations(TCSlice &slc, const CTP_t &inCTP)
T get(std::string const &key) const
std::vector< float > neutralVxCuts
void MakeHaloTj(TCSlice &slc, Trajectory &muTj, bool prt)
std::vector< short > EndVtx
std::vector< VtxStore > vtxs
2D vertices
std::array< unsigned short, 2 > EndPt
First and last point in the trajectory that has charge.
tick_as<> tick
Tick number, represented by std::ptrdiff_t.
std::vector< float > match3DCuts
3D matching cuts
void ScoreVertices(TCSlice &slc)
std::string PrintHit(const TCHit &tch)
std::vector< float > Envelope
call study functions to develop cuts, etc
const geo::GeometryCore * geom
std::vector< unsigned short > minPts
min number of Pts required to make a trajectory
bool has_key(std::string const &key) const
PlaneID_t Plane
Index of the plane within its TPC.
void ReleaseHits(TCSlice &slc, Trajectory const &tj)
void DefineTjParents(TCSlice &slc, bool prt)
float HitsPosTick(const TCSlice &slc, const std::vector< unsigned int > &hitsInMultiplet, float &sum, HitStatus_t hitRequest)
std::vector< unsigned int > firstWire
the first wire with a hit
std::vector< float > BeginChg
int ID
ID that is local to one slice.
std::vector< TCSlice > slices
Detector simulation of raw signals on wires.
std::bitset< 16 > modes
number of points to find AveChg
std::vector< TCHit > slHits
std::vector< float > chargeCuts
bool aveHitRMSValid
set true when the average hit RMS is well-known
cet::coded_exception< errors::ErrorCodes, ExceptionDetail::translate > Exception
std::vector< int > EnvPlane
std::vector< TrajPoint > seeds
float maxWireSkipNoSignal
max number of wires to skip w/o a signal on them
std::vector< short > MCSMom
void CheckTraj(TCSlice &slc, Trajectory &tj)
std::vector< float > vtxScoreWeights
decltype(auto) get(T &&obj)
ADL-aware version of std::to_string.
bool DecodeDebugString(std::string strng)
std::vector< int > StageNum
std::vector< recob::Hit > const * srcHits
std::bitset< 128 > useAlg
Max hit separation for making junk trajectories. < 0 to turn off.
Contains all timing reference information for the detector.
std::vector< float > pfpStitchCuts
cuts for stitching between TPCs
std::vector< short > deltaRayTag
short StepDir
-1 = going US (-> small wire#), 1 = going DS (-> large wire#)
std::bitset< 128 > AlgMod
Bit set if algorithm AlgBit_t modifed the trajectory.
call study functions to develop cuts, etc
void LastEndMerge(TCSlice &slc, CTP_t inCTP)
std::vector< short > muonTag
void UpdateVxEnvironment(TCSlice &slc)
std::vector< float > BeginWir
geo::PlaneID DecodeCTP(CTP_t CTP)
std::optional< T > get_if_present(std::string const &key) const
std::vector< float > EndChg
TrajClusterAlg(fhicl::ParameterSet const &pset)
std::bitset< 128 > dbgAlg
Allow user to turn on debug printing in algorithms (that print...)
std::vector< recob::Hit > const * allHits
CTP_t EncodeCTP(unsigned int cryo, unsigned int tpc, unsigned int plane)
std::vector< unsigned int > nWires
void KillPoorVertices(TCSlice &slc)
std::vector< float > chkStopCuts
Bragg peak finder configuration.
int ID
ID of the recob::Slice (not the sub-slice)
void EndMerge(TCSlice &slc, CTP_t inCTP, bool lastPass)
void MergeTPHits(std::vector< unsigned int > &tpHits, std::vector< recob::Hit > &newHitCol, std::vector< unsigned int > &newHitAssns) const
unsigned int Nplanes(TPCID const &tpcid=tpc_zero) const
Returns the total number of planes in the specified TPC.
std::vector< float > vtx2DCuts
Max position pull, max Position error rms.
std::bitset< 8 > Strategy
void Finish3DShowers(TCSlice &slc)
IDparameter< geo::TPCID > TPCID
Member type of validated geo::TPCID parameter.
void SetSourceHits(std::vector< recob::Hit > const &srcHits)
void FindPFParticles(detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, TCSlice &slc)
float multHitSep
preferentially "merge" hits with < this separation
2D representation of charge deposited in the TDC/wire plane
void ChkInTraj(std::string someText, TCSlice &slc)
void ReconstructAllTraj(detinfo::DetectorPropertiesData const &detProp, TCSlice &slc, CTP_t inCTP)
unsigned short NumPtsWithCharge(const TCSlice &slc, const Trajectory &tj, bool includeDeadWires)
std::vector< int > EnvShowerID
void MakeJunkVertices(TCSlice &slc, const CTP_t &inCTP)
call study functions to develop cuts, etc (see TCTruth.cxx)
void Find2DVertices(detinfo::DetectorPropertiesData const &detProp, TCSlice &slc, const CTP_t &inCTP, unsigned short pass)
std::vector< int > IsShowerTj
std::vector< short > PlaneNum
short recoSlice
only reconstruct the slice with ID (0 = all)
master switch for turning on debug mode
void DefinePFPParents(TCSlice &slc)
don't mess with this line
art framework interface to geometry description
CTP_t CTP
set to an invalid CTP
Event finding and building.
void SetTPEnvironment(TCSlice &slc, CTP_t inCTP)
void FindJunkTraj(TCSlice &slc, CTP_t inCTP)