37 : fCaloAlg(pset.
get<
fhicl::ParameterSet>(
"CaloAlg")), fMVAReader(
"Silent")
41 bool badinput =
false;
48 if (pset.
has_key(
"Mode")) userMode = pset.
get<
short>(
"Mode");
52 if (pset.
has_key(
"StudyMode")) {
53 std::cout <<
"StudyMode is not valid anymore. Try Study: 1 or Study: 2, etc/n";
56 userMode = pset.
get<
short>(
"Study");
62 if (pset.
has_key(
"SaveShowerTree"))
66 std::vector<std::string> skipAlgsVec;
67 if (pset.
has_key(
"SkipAlgs")) skipAlgsVec = pset.
get<std::vector<std::string>>(
"SkipAlgs");
68 std::vector<std::string> debugConfigVec;
69 if (pset.
has_key(
"DebugConfig"))
70 debugConfigVec = pset.
get<std::vector<std::string>>(
"DebugConfig");
74 std::vector<float> aveHitRMS;
75 if (pset.
has_key(
"AveHitRMS")) aveHitRMS = pset.
get<std::vector<float>>(
"AveHitRMS");
77 if (!aveHitRMS.empty()) {
84 tcc.
minPts = pset.
get<std::vector<unsigned short>>(
"MinPts");
88 tcc.
chargeCuts = pset.
get<std::vector<float>>(
"ChargeCuts", {3, 0.15, 0.25});
98 tcc.
muonTag = pset.
get<std::vector<short>>(
"MuonTag", {-1, -1, -1, -1});
100 tcc.
showerTag = pset.
get<std::vector<float>>(
"ShowerTag", {-1, -1, -1, -1, -1, -1});
101 std::string fMVAShowerParentWeights =
"NA";
102 pset.
get_if_present<std::string>(
"MVAShowerParentWeights", fMVAShowerParentWeights);
104 if (pset.
has_key(
"MatchTruth")) {
105 std::cout <<
"MatchTruth is not used. Use ClusterAnaV2 or DebugConfig to configure\n";
107 tcc.
vtx2DCuts = pset.
get<std::vector<float>>(
"Vertex2DCuts", {-1, -1, -1, -1, -1, -1, -1});
108 tcc.
vtx3DCuts = pset.
get<std::vector<float>>(
"Vertex3DCuts", {-1, -1});
110 tcc.
match3DCuts = pset.
get<std::vector<float>>(
"Match3DCuts", {-1, -1, -1, -1, -1});
123 <<
"Bad input from fcl file. Vector lengths for MinPtsFit, MaxAngleRange and MinMCSMom " 124 "should be defined for each reconstruction pass";
128 <<
"Vertex2DCuts must be size 11\n 0 = Max length definition for short TJs\n 1 = Max " 129 "vtx-TJ sep short TJs\n 2 = Max vtx-TJ sep long TJs\n 3 = Max position pull for >2 " 130 "TJs\n 4 = Max vtx position error\n 5 = Min MCSMom for one of two TJs\n 6 = Min " 131 "fraction of wires hit btw vtx and Tjs\n 7 = Min Score\n 8 = min ChgFrac at a vtx or " 132 "merge point\n 9 = max MCSMom asymmetry, 10 = require chg on wires btw vtx and tjs in " 141 <<
"Vertex3DCuts must be size > 2\n 0 = 2D Vtx max dX (cm)\n 1 = 2D Vtx max pull\n 2 = max " 142 "3D separation (cm) btw PFP and vertex";
144 std::cout <<
"WARNING: Vertex3DCuts is size 2 but should be size 3, where Vertex3DCuts[2] = " 145 "max 3D separation (cm) btw a PFP and a 3D vertex. Setting it to 3 cm\n";
150 <<
"KinkCuts must be size 3\n 0 = Number of points to fit at the end of the trajectory\n 1 " 151 "= Minimum kink significance\n 2 = Use charge in significance calculation? (yes if > " 152 "0)\n 3 = 3D kink fit length (cm). \nYou are using an out-of-date specification?\n";
158 <<
"Are you using an out-of-date specification for KinkCuts? KinkCuts[0] is the number of " 164 <<
"ChargeCuts must be size 3\n 0 = Charge pull cut\n 1 = Min allowed fractional chg RMS\n " 165 "2 = Max allowed fractional chg RMS";
169 <<
"MuonTag must be size 4\n 0 = minPtsFit\n 1 = minMCSMom\n 2= maxWireSkipNoSignal\n 3 = " 170 "min delta ray length for tagging\n 4 = dress muon window size (optional)";
173 <<
"DeltaRayTag must be size 3\n 0 = Max endpoint sep\n 1 = min MCSMom\n 2 = max MCSMom";
176 <<
"ChkStopCuts must be size 3\n 0 = Min Charge ratio\n 1 = Charge slope pull cut\n 2 = " 177 "Charge fit chisq cut\n 3 = BraggSplit FOM (optional)";
180 <<
"ShowerTag must be size 13\n 0 = Mode\n 1 = max MCSMom\n 2 = max separation (WSE " 181 "units)\n 3 = Max angle diff\n 4 = Factor * rms width\n 5 = Min half width\n 6 = min " 182 "total Tps\n 7 = Min Tjs\n 8 = max parent FOM\n 9 = max direction FOM 10 = max " 183 "AspectRatio\n 11 = min Score to preserve a vertex\n 12 = Debug showers in CTP\n";
184 std::cout <<
" Fixing this problem...";
193 <<
"Match3DCuts must be size 5\n 0 = dx (cm) matching cut\n 1 = max number of 3D " 194 "combinations\n 2 = min length for 2-view match\n 3 = number of TP3Ds in each plane to " 195 "fit in each PFP section\n 4 = max pull for accepting TP3Ds in sections\n 5 = max " 196 "ChiDOF for a SectionFit";
200 mf::LogVerbatim(
"TC") <<
"Last element of AngleRange != 90 degrees. Fixing it\n";
216 for (
auto strng : debugConfigVec) {
220 for (
unsigned short ib = 0; ib <
AlgBitNames.size(); ++ib) {
229 std::cout <<
"DecodeDebugString failed: " << strng <<
"\n";
237 if (range < 0 || range > 90)
239 <<
"Invalid angle range " << range <<
" Must be 0 - 90 degrees";
249 <<
"Increase the size of UseAlgs to at least " <<
kAlgBitSize;
258 <<
"Increase the size of EndFlag to at least " <<
kFlagBitSize;
260 bool printHelp =
false;
261 for (
unsigned short ib = 0; ib <
AlgBitNames.size(); ++ib)
268 for (
auto strng : skipAlgsVec) {
270 if (strng ==
"All") {
272 for (
unsigned short ib = 0; ib <
AlgBitNames.size(); ++ib)
277 for (
unsigned short ib = 0; ib <
AlgBitNames.size(); ++ib) {
285 std::cout <<
"******* Unknown SkipAlgs input string '" << strng <<
"'\n";
290 std::cout <<
"Valid AlgNames:";
292 std::cout <<
" " << strng;
294 std::cout <<
"Or specify All to turn all algs off\n";
297 if (fMVAShowerParentWeights !=
"NA" &&
tcc.
showerTag[0] > 0)
317 tcc.
geom = lar::providerFrom<geo::Geometry>();
340 thr = {UINT_MAX, UINT_MAX};
343 unsigned int tpc =
hit.WireID().TPC;
353 std::vector<unsigned int>& hitsInSlice,
359 if (hitsInSlice.size() < 2)
return;
362 if (!
CreateSlice(clockData, detProp, hitsInSlice, sliceID))
return;
375 <<
" AveHitRMS vector size != the number of planes ";
377 std::cout <<
"Reconstruct " << hitsInSlice.size() <<
" hits in Slice " << sliceID
378 <<
" in TPC " << slc.TPCID.TPC <<
"\n";
379 for (
unsigned short plane = 0; plane < slc.nPlanes; ++plane) {
380 CTP_t inCTP =
EncodeCTP(slc.TPCID.Cryostat, slc.TPCID.TPC, plane);
382 if (!slc.isValid)
return;
391 for (
unsigned short plane = 0; plane < slc.nPlanes; ++plane) {
392 CTP_t inCTP =
EncodeCTP(slc.TPCID.Cryostat, slc.TPCID.TPC, plane);
394 std::cout <<
"RTC: ChkVtxAssociations found an error\n";
407 std::cout <<
"SHOWER TREE STAGE NUM SIZE: " <<
stv.
StageNum.size() << std::endl;
413 mf::LogVerbatim(
"TC") <<
"RunTrajCluster failed in MakeAllTrajClusters";
423 for (
auto& tj : slc.tjs) {
424 for (
unsigned short ib = 0; ib <
AlgBitNames.size(); ++ib)
429 slc.mallTraj.resize(0);
443 if (nwires > slc.
nWires[plane])
return;
449 for (
unsigned short pass = 0; pass <
tcc.
minPtsFit.size(); ++pass) {
450 for (
unsigned int ii = 0; ii < nwires; ++ii) {
452 unsigned int iwire = 0;
453 unsigned int jwire = 0;
461 iwire = slc.
lastWire[plane] - ii - 1;
464 if (iwire > slc.
wireHitRange[plane].size() - 1)
continue;
465 if (jwire > slc.
wireHitRange[plane].size() - 1)
continue;
467 if (slc.
wireHitRange[plane][iwire].first == UINT_MAX)
continue;
468 if (slc.
wireHitRange[plane][jwire].first == UINT_MAX)
continue;
469 unsigned int ifirsthit = slc.
wireHitRange[plane][iwire].first;
470 unsigned int ilasthit = slc.
wireHitRange[plane][iwire].second;
471 unsigned int jfirsthit = slc.
wireHitRange[plane][jwire].first;
472 unsigned int jlasthit = slc.
wireHitRange[plane][jwire].second;
473 if (ifirsthit > slc.
slHits.size() || ilasthit > slc.
slHits.size()) {
474 std::cout <<
"RAT: bad hit range " << ifirsthit <<
" " << ilasthit <<
" size " 475 << slc.
slHits.size() <<
" inCTP " << inCTP <<
"\n";
478 for (
unsigned int iht = ifirsthit; iht <= ilasthit; ++iht) {
481 mf::LogVerbatim(
"TC") <<
"+++++++ Pass " << pass <<
" Found debug hit " 486 if (slc.
slHits[iht].InTraj != 0)
continue;
490 std::cout <<
"RAT: Bad allHitsIndex\n";
495 unsigned int fromWire = iHit.WireID().Wire;
496 float fromTick = iHit.PeakTime();
497 float iqtot = iHit.Integral();
498 float hitsRMSTick = iHit.RMS();
499 std::vector<unsigned int> iHitsInMultiplet;
502 iHitsInMultiplet.resize(1);
503 iHitsInMultiplet[0] = iht;
507 if (iHitsInMultiplet.empty())
continue;
509 if (iHitsInMultiplet.size() > 4)
continue;
511 if (iHitsInMultiplet.size() > 1) {
515 if (hitsRMSTick == 0)
continue;
516 bool fatIHit = (hitsRMSTick > maxHitsRMS);
518 mf::LogVerbatim(
"TC") <<
" hit RMS " << iHit.RMS() <<
" BB Multiplicity " 519 << iHitsInMultiplet.size() <<
" AveHitRMS[" << plane <<
"] " 520 <<
evt.
aveHitRMS[plane] <<
" HitsRMSTick " << hitsRMSTick
521 <<
" fatIHit " << fatIHit;
522 for (
unsigned int jht = jfirsthit; jht <= jlasthit; ++jht) {
524 if (slc.
slHits[iht].InTraj != 0)
break;
525 if (slc.
slHits[jht].InTraj != 0)
continue;
527 for (
auto& slHit : slc.
slHits) {
528 if (slHit.InTraj < 0) {
529 std::cout <<
"RAT: Dirty hit " <<
PrintHit(slHit) <<
" EventsProcessed " 534 unsigned int toWire = jwire;
537 float toTick = jHit.PeakTime();
538 float jqtot = jHit.Integral();
539 std::vector<unsigned int> jHitsInMultiplet;
542 jHitsInMultiplet.resize(1);
543 jHitsInMultiplet[0] = jht;
547 if (jHitsInMultiplet.empty())
continue;
549 if (jHitsInMultiplet.size() > 4)
continue;
552 if (hitsRMSTick == 0)
continue;
553 bool fatJHit = (hitsRMSTick > maxHitsRMS);
556 if ((fatIHit && !fatJHit) || (!fatIHit && fatJHit)) {
continue; }
560 if (jHitsInMultiplet.size() > 1)
563 bool hitsOK =
TrajHitsOK(slc, iHitsInMultiplet, jHitsInMultiplet);
565 if (!hitsOK)
continue;
568 if (!
StartTraj(work, fromWire, fromTick, toWire, toTick, inCTP, pass))
continue;
571 std::cout <<
"RAT: StartTraj major failure\n";
574 if (work.
Pts.empty()) {
587 std::cout <<
"RAT: AddHits major failure\n";
590 if (!sigOK || work.
Pts[0].Chg == 0) {
596 work.
Pts[0].Pos = work.
Pts[0].HitPos;
605 std::cout <<
"RAT: StepAway major failure\n";
614 std::cout <<
"RAT: CheckTraj major failure\n";
619 <<
"-" << work.
EndPt[1] <<
" IsGood " << work.
IsGood;
622 <<
"StepAway done: IsGood " << work.
IsGood <<
" NumPtsWithCharge " 635 auto& tj = slc.
tjs[slc.
tjs.size() - 1];
638 std::cout <<
"RAT: InTrajOK major failure. " << tj.ID <<
"\n";
651 for (
auto tp :
seeds) {
652 unsigned short nAvailable = 0;
653 for (
unsigned short ii = 0; ii < tp.Hits.size(); ++ii) {
654 if (!tp.UseHit[ii])
continue;
655 unsigned int iht = tp.Hits[ii];
656 if (slc.
slHits[iht].InTraj == 0) ++nAvailable;
663 if (nAvailable == 0)
continue;
666 for (
unsigned short ii = 0; ii < tp.Hits.size(); ++ii) {
667 if (!tp.UseHit[ii])
continue;
668 unsigned int iht = tp.Hits[ii];
669 if (slc.
slHits[iht].InTraj == 0) slc.
slHits[iht].InTraj = work.
ID;
673 if (tp.Dir[0] < 0) work.
StepDir = -1;
680 work.
Pts.push_back(tp);
685 std::cout <<
"RAT: CheckTraj major failure\n";
699 auto& tj = slc.
tjs[slc.
tjs.size() - 1];
724 for (
auto& tj : slc.
tjs) {
725 if (tj.AlgMod[
kKilled])
continue;
726 if (tj.PDGCode != 13)
continue;
741 std::cout <<
"RAT: MakeJunkVertices major failure\n";
746 for (
unsigned short ivx = 0; ivx < slc.
vtxs.size(); ++ivx) {
747 auto& vx2 = slc.
vtxs[ivx];
748 if (vx2.ID == 0)
continue;
749 if (vx2.CTP != inCTP)
continue;
761 std::cout <<
"RAT: ChkVtxAssociations found an error. Events processed " 778 for (
auto& slHit : slc.
slHits)
779 if (slHit.InTraj < 0) slHit.InTraj = 0;
783 std::vector<unsigned int> tHits;
785 for (
unsigned int iwire = slc.
firstWire[plane]; iwire < slc.
lastWire[plane] - 3; ++iwire) {
789 unsigned int jwire = iwire + 1;
792 unsigned int ifirsthit = slc.
wireHitRange[plane][iwire].first;
793 unsigned int ilasthit = slc.
wireHitRange[plane][iwire].second;
794 unsigned int jfirsthit = slc.
wireHitRange[plane][jwire].first;
795 unsigned int jlasthit = slc.
wireHitRange[plane][jwire].second;
796 for (
unsigned int iht = ifirsthit; iht <= ilasthit; ++iht) {
797 if (iht >= slc.
slHits.size())
break;
798 auto& islHit = slc.
slHits[iht];
799 if (islHit.InTraj != 0)
continue;
800 std::vector<unsigned int> iHits;
804 if (prt)
mf::LogVerbatim(
"TC") <<
"FJT: debug iht multiplet size " << iHits.size();
805 if (iHits.empty())
continue;
806 for (
unsigned int jht = jfirsthit; jht <= jlasthit; ++jht) {
807 if (jht >= slc.
slHits.size())
break;
808 auto& jslHit = slc.
slHits[jht];
809 if (jslHit.InTraj != 0)
continue;
810 if (prt &&
HitSep2(slc, iht, jht) < 100)
814 std::vector<unsigned int> jHits;
816 if (jHits.empty())
continue;
821 for (
auto iht : iHits)
822 if (slc.
slHits[iht].InTraj == 0) tHits.push_back(iht);
823 for (
auto jht : jHits)
824 if (slc.
slHits[jht].InTraj == 0) tHits.push_back(jht);
825 for (
auto tht : tHits)
826 slc.
slHits[tht].InTraj = -4;
828 if (iwire != 0) { loWire = iwire - 1; }
832 unsigned int hiWire = jwire + 1;
833 if (hiWire > slc.
nWires[plane])
break;
834 unsigned short nit = 0;
836 bool hitsAdded =
false;
837 for (
unsigned int kwire = loWire; kwire <= hiWire; ++kwire) {
838 if (slc.
wireHitRange[plane][kwire].first == UINT_MAX)
continue;
839 if (slc.
wireHitRange[plane][kwire].second == UINT_MAX)
continue;
840 unsigned int kfirsthit = slc.
wireHitRange[plane][kwire].first;
841 unsigned int klasthit = slc.
wireHitRange[plane][kwire].second;
842 for (
unsigned int kht = kfirsthit; kht <= klasthit; ++kht) {
843 if (kht >= slc.
slHits.size())
continue;
844 if (slc.
slHits[kht].InTraj != 0)
continue;
846 if (std::find(tHits.begin(), tHits.end(), kht) != tHits.end())
continue;
851 for (
auto jht : jHits) {
852 if (slc.
slHits[jht].InTraj != 0)
continue;
853 tHits.push_back(jht);
854 slc.
slHits[jht].InTraj = -4;
855 if (kwire > hiWire) hiWire = kwire;
856 if (kwire < loWire) loWire = kwire;
861 if (!hitsAdded)
break;
864 if (hiWire >= slc.
nWires[plane])
break;
867 for (
auto iht : tHits)
868 slc.
slHits[iht].InTraj = 0;
869 if (tHits.size() < 3)
continue;
872 myprt <<
"FJT: tHits";
873 for (
auto tht : tHits)
878 if (
IsGhost(slc, tHits))
break;
883 if (slc.
slHits[jht].InTraj > 0)
break;
900 unsigned short itj = 0;
901 std::vector<unsigned int> tHits;
902 std::vector<unsigned int> atHits;
903 for (
auto& tj : slc.
tjs) {
905 if (tj.AlgMod[
kKilled])
continue;
907 for (
auto& tp : tj.Pts) {
908 if (tp.Hits.size() > 16) {
911 <<
"ChkInTraj: More than 16 hits created a UseHit bitset overflow\n";
913 std::cout <<
"ChkInTraj major failure\n";
918 std::cout << someText <<
" ChkInTraj hit size mis-match in tj ID " << tj.ID
920 for (
unsigned short ib = 0; ib <
AlgBitNames.size(); ++ib)
921 if (tj.AlgMod[ib]) std::cout <<
" " <<
AlgBitNames[ib];
926 if (tHits.size() < 2) {
927 mf::LogVerbatim(
"TC") << someText <<
" ChkInTraj: Insufficient hits in traj " << tj.ID;
932 std::sort(tHits.begin(), tHits.end());
934 for (iht = 0; iht < slc.
slHits.size(); ++iht) {
935 if (slc.
slHits[iht].InTraj == tID) atHits.push_back(iht);
937 if (atHits.size() < 2) {
938 mf::LogVerbatim(
"TC") << someText <<
" ChkInTraj: Insufficient hits in atHits in traj " 939 << tj.ID <<
" Killing it";
943 if (!std::equal(tHits.begin(), tHits.end(), atHits.begin())) {
945 myprt << someText <<
" ChkInTraj failed: inTraj - UseHit mis-match for tj ID " << tID
946 <<
" tj.WorkID " << tj.WorkID <<
" atHits size " << atHits.size() <<
" tHits size " 947 << tHits.size() <<
" in CTP " << tj.CTP <<
"\n";
948 myprt <<
"AlgMods: ";
949 for (
unsigned short ib = 0; ib <
AlgBitNames.size(); ++ib)
950 if (tj.AlgMod[ib]) myprt <<
" " <<
AlgBitNames[ib];
952 myprt <<
"index inTraj UseHit \n";
953 for (iht = 0; iht < atHits.size(); ++iht) {
954 myprt <<
"iht " << iht <<
" " <<
PrintHit(slc.
slHits[atHits[iht]]);
955 if (iht < tHits.size()) myprt <<
" " <<
PrintHit(slc.
slHits[tHits[iht]]);
956 if (atHits[iht] != tHits[iht]) myprt <<
" <<< " << atHits[iht] <<
" != " << tHits[iht];
960 if (tHits.size() > atHits.size()) {
961 for (iht = atHits.size(); iht < atHits.size(); ++iht) {
962 myprt <<
"atHits " << iht <<
" " <<
PrintHit(slc.
slHits[atHits[iht]]) <<
"\n";
968 for (
unsigned short end = 0;
end < 2; ++
end) {
969 if (tj.VtxID[
end] > slc.
vtxs.size()) {
970 mf::LogVerbatim(
"TC") << someText <<
" ChkInTraj: Bad VtxID " << tj.ID;
971 std::cout << someText <<
" ChkInTraj: Bad VtxID " << tj.ID <<
" vtx size " 972 << slc.
vtxs.size() <<
"\n";
987 std::vector<recob::Hit>& newHitCol,
988 std::vector<unsigned int>& newHitAssns)
const 993 if (tpHits.empty())
return;
996 if (tpHits.size() == 1) {
998 std::cout <<
"MergeTPHits Bad input hit index " << tpHits[0] <<
" allHits size " 1003 newHitAssns[tpHits[0]] = newHitCol.size() - 1;
1008 std::vector<unsigned int> wires;
1009 std::vector<std::vector<unsigned int>> wireHits;
1011 wires.push_back(firstHit.WireID().Wire);
1012 std::vector<unsigned int>
tmp(1, tpHits[0]);
1013 wireHits.push_back(tmp);
1014 for (
unsigned short ii = 1; ii < tpHits.size(); ++ii) {
1016 unsigned int wire =
hit.WireID().Wire;
1017 unsigned short indx = 0;
1018 for (indx = 0; indx < wires.size(); ++indx)
1019 if (wires[indx] == wire)
break;
1020 if (indx == wires.size()) {
1021 wires.push_back(wire);
1022 wireHits.resize(wireHits.size() + 1);
1024 wireHits[indx].push_back(tpHits[ii]);
1028 for (
unsigned short indx = 0; indx < wireHits.size(); ++indx) {
1029 auto& hitsOnWire = wireHits[indx];
1031 for (
unsigned short ii = 0; ii < hitsOnWire.size(); ++ii) {
1032 newHitAssns[hitsOnWire[ii]] = newHitCol.size() - 1;
1048 if (tpHits.size() == 1) {
1050 std::cout <<
"MergeTPHits Bad input hit index " << tpHits[0] <<
" allHits size " 1062 oldHit.SigmaPeakTime(),
1064 oldHit.PeakAmplitude(),
1065 oldHit.SigmaPeakAmplitude(),
1066 oldHit.ROISummedADC(),
1067 oldHit.HitSummedADC(),
1069 oldHit.SigmaIntegral(),
1075 oldHit.SignalType(),
1079 double integral = 0;
1080 double sIntegral = 0;
1081 double peakTime = 0;
1082 double sPeakTime = 0;
1084 double sPeakAmp = 0;
1085 float ROIsumADC = 0;
1086 float HitsumADC = 0;
1089 for (
auto allHitsIndex : tpHits) {
1092 if (
hit.StartTick() < startTick) startTick =
hit.StartTick();
1093 if (
hit.EndTick() > endTick) endTick =
hit.EndTick();
1094 double intgrl =
hit.Integral();
1095 sPeakTime += intgrl *
hit.SigmaPeakTime();
1096 sPeakAmp += intgrl *
hit.SigmaPeakAmplitude();
1097 ROIsumADC +=
hit.ROISummedADC();
1098 HitsumADC +=
hit.HitSummedADC();
1100 sIntegral += intgrl *
hit.SigmaIntegral();
1103 if (integral <= 0) {
1104 std::cout <<
"MergeTPHits found bad hit integral " << integral <<
"\n";
1109 std::vector<double> shape(endTick - startTick + 1, 0.);
1110 for (
auto allHitsIndex : tpHits) {
1112 double peakTick =
hit.PeakTime();
1113 double rms =
hit.RMS();
1114 double peakAmp =
hit.PeakAmplitude();
1116 double arg = ((double)
tick - peakTick) / rms;
1117 unsigned short indx =
tick - startTick;
1118 shape[indx] += peakAmp * exp(-0.5 * arg * arg);
1126 unsigned short indx =
tick - startTick;
1127 sigsum += shape[indx];
1128 sigsumt += shape[indx] *
tick;
1131 peakTime = sigsumt / sigsum;
1133 sPeakTime /= integral;
1138 double dTick =
tick - peakTime;
1139 unsigned short indx =
tick - startTick;
1140 sigsumt += shape[indx] * dTick * dTick;
1142 double rms = std::sqrt(sigsumt / sigsum);
1145 double chgNorm = 2.507 * firstHit.PeakAmplitude() * firstHit.RMS() / firstHit.Integral();
1147 peakAmp = (float)(integral * chgNorm / (2.507 * rms));
1149 sPeakAmp /= integral;
1151 startTick = peakTime - 3 * rms;
1152 endTick = peakTime + 3 * rms;
1172 firstHit.SignalType(),
1221 std::vector<unsigned int>& hitsInSlice,
1227 if (hitsInSlice.size() < 2)
return false;
1231 slc.
slHits.resize(hitsInSlice.size());
1233 unsigned int cstat = 0;
1234 unsigned int tpc = UINT_MAX;
1235 unsigned int cnt = 0;
1236 std::vector<unsigned int> nHitsInPln;
1237 for (
auto iht : hitsInSlice) {
1238 if (iht >= (*
evt.
allHits).size())
return false;
1241 cstat =
hit.WireID().Cryostat;
1242 tpc =
hit.WireID().TPC;
1247 if (
hit.WireID().Cryostat != cstat ||
hit.WireID().TPC != tpc)
return false;
1248 slc.
slHits[cnt].allHitsIndex = iht;
1249 ++nHitsInPln[
hit.WireID().Plane];
1253 for (
auto hip : nHitsInPln)
1254 if (hip < 2)
return false;
1264 if (
tcc.
dbgSlc) std::cout <<
"Enabled debugging in sub-slice " <<
slices.size() - 1 <<
"\n";
1272 for (
unsigned int iht = 0; iht < slc.
slHits.size(); ++iht) {
1273 unsigned int ahi = slc.
slHits[iht].allHitsIndex;
1275 std::cout <<
"CreateSlice: Bad allHitsIndex " << ahi <<
" " << (*
evt.
allHits).
size()
1290 for (
auto& slc :
slices) {
1291 if (!slc.isValid)
continue;
1293 for (
auto& pfp : slc.pfps) {
1294 if (pfp.ID <= 0)
continue;
1295 pfp.TjUIDs.resize(pfp.TjIDs.size());
1296 for (
unsigned short ii = 0; ii < pfp.TjIDs.size(); ++ii) {
1298 if (pfp.TjIDs[ii] <= 0 || pfp.TjIDs[ii] > (
int)slc.tjs.size()) {
1299 std::cout <<
"FinishEvent found an invalid T" << pfp.TjIDs[ii] <<
" in P" << pfp.UID
1301 slc.isValid =
false;
1304 auto& tj = slc.tjs[pfp.TjIDs[ii] - 1];
1305 pfp.TjUIDs[ii] = tj.UID;
1312 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)
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.
const geo::WireReadoutGeom * wireReadoutGeom
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)
unsigned int NTPC(CryostatID const &cryoid=details::cryostat_zero) const
Returns the total number of TPCs in the specified cryostat.
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)
cout<< "Opened file "<< fin<< " ixs= "<< ixs<< endl;if(ixs==0) hhh=(TH1F *) fff-> Get("h1")
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.
unsigned int Nplanes(TPCID const &tpcid=details::tpc_zero) const
Returns the total number of planes in the specified TPC.
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
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)