54 std::vector<std::vector<unsigned int>> matchedclusters;
61 std::vector<std::vector<TH1D>> signals(nplanes);
63 std::vector<std::vector<unsigned int>> Cls(nplanes);
64 std::vector<std::vector<CluLen>> clulens(nplanes);
66 for (
size_t iclu = 0; iclu < clusterlist.size(); ++iclu) {
70 float w0 = clusterlist[iclu]->StartWire();
71 float w1 = clusterlist[iclu]->EndWire();
72 float t0 = clusterlist[iclu]->StartTick();
73 float t1 = clusterlist[iclu]->EndTick();
80 clusterlist[iclu]->
Plane().TPC,
81 clusterlist[iclu]->
Plane().Cryostat);
84 clusterlist[iclu]->
Plane().TPC,
85 clusterlist[iclu]->
Plane().Cryostat);
86 clulen.length = std::hypot((w0 - w1) * wire_pitch, x0 -
x1);
88 switch (clusterlist[iclu]->View()) {
90 if (
fEnableU) clulens[0].push_back(clulen);
93 if (
fEnableV) clulens[1].push_back(clulen);
96 if (
fEnableZ) clulens[2].push_back(clulen);
103 for (
size_t i = 0; i < clulens.size(); ++i) {
104 std::sort(clulens[i].
begin(), clulens[i].
end(), SortByLength);
105 for (
size_t j = 0; j < clulens[i].size(); ++j) {
106 Cls[i].push_back(clulens[i][j].index);
110 for (
int i = 0; i < nplanes; ++i) {
111 for (
size_t ic = 0; ic < Cls[i].size(); ++ic) {
113 Form(
"sig_%d_%d", i,
int(ic)), Form(
"sig_%d_%d", i,
int(ic)), nts + 100, -100, nts);
115 Form(
"sigint_%d_%d", i,
int(ic)), Form(
"sigint_%d_%d", i,
int(ic)), nts + 100, -100, nts);
116 std::vector<art::Ptr<recob::Hit>> hitlist = fm.at(Cls[i][ic]);
117 std::sort(hitlist.begin(), hitlist.end(),
SortByWire);
119 for (
auto theHit = hitlist.begin(); theHit != hitlist.end(); theHit++) {
121 double time = (*theHit)->PeakTime();
123 (*theHit)->WireID().Plane, (*theHit)->WireID().TPC, (*theHit)->WireID().Cryostat);
125 double charge = (*theHit)->Integral();
126 int bin = sig.FindBin(time);
127 sig.SetBinContent(bin, sig.GetBinContent(bin) + charge);
128 for (
int j = bin; j <= sig.GetNbinsX(); ++j) {
129 sigint.SetBinContent(j, sigint.GetBinContent(j) + charge);
132 if (sigint.Integral()) sigint.Scale(1. / sigint.GetBinContent(sigint.GetNbinsX()));
133 signals[i].push_back(sigint);
138 std::vector<int> matched(clusterlist.size());
139 for (
size_t i = 0; i < clusterlist.size(); ++i)
142 for (
int i = 0; i < nplanes - 1; ++i) {
143 for (
int j = i + 1; j < nplanes; ++j) {
144 for (
size_t c1 = 0;
c1 < Cls[i].size(); ++
c1) {
145 for (
size_t c2 = 0;
c2 < Cls[j].size(); ++
c2) {
148 if (clusterlist[Cls[i][
c1]]->View() == clusterlist[Cls[j][
c2]]->View())
continue;
150 if (clusterlist[Cls[i][
c1]]->
Plane().Cryostat !=
151 clusterlist[Cls[j][
c2]]->Plane().Cryostat)
153 if (clusterlist[Cls[i][
c1]]->
Plane().TPC != clusterlist[Cls[j][
c2]]->Plane().TPC)
156 if (matched[Cls[i][
c1]] == 1 && matched[Cls[j][
c2]] == 1)
continue;
159 if (signals[i][
c1].Integral() && signals[j][
c2].Integral())
160 ks = signals[i][
c1].KolmogorovTest(&signals[j][
c2]);
163 <<
"One of the two clusters appears to be empty: " << clusterlist[Cls[i][
c1]]->ID()
164 <<
" " << clusterlist[Cls[j][
c2]]->ID() <<
" " << i <<
" " << j <<
" " <<
c1 <<
" " 165 << c2 <<
" " << signals[i][
c1].Integral() <<
" " << signals[j][
c2].Integral();
174 for (
size_t l = 0; l < matchedclusters.size(); ++l) {
175 for (
size_t m = 0; m < matchedclusters[l].size(); ++m) {
176 if (matchedclusters[l][m] == Cls[i][
c1]) {
180 else if (matchedclusters[l][m] == Cls[j][c2]) {
188 bool matchview =
false;
190 for (
size_t ii = 0; ii < matchedclusters[imatch].size(); ++ii) {
191 if (clusterlist[matchedclusters[imatch][ii]]->View() ==
192 clusterlist[Cls[i][
c1]]->View()) {
195 if (fm.at(Cls[i][
c1]).size() > fm.at(matchedclusters[imatch][ii]).size()) {
196 matched[matchedclusters[imatch][ii]] = 0;
197 matchedclusters[imatch][ii] = Cls[i][
c1];
198 matched[Cls[i][
c1]] = 1;
203 matchedclusters[imatch].push_back(Cls[i][c1]);
204 matched[Cls[i][
c1]] = 1;
208 bool matchview =
false;
209 for (
size_t jj = 0; jj < matchedclusters[imatch].size(); ++jj) {
210 if (clusterlist[matchedclusters[imatch][jj]]->View() ==
211 clusterlist[Cls[j][
c2]]->View()) {
214 if (fm.at(Cls[j][c2]).size() > fm.at(matchedclusters[imatch][jj]).size()) {
215 matched[matchedclusters[imatch][jj]] = 0;
216 matchedclusters[imatch][jj] = Cls[j][
c2];
217 matched[Cls[j][
c2]] = 1;
222 matchedclusters[imatch].push_back(Cls[j][c2]);
223 matched[Cls[j][
c2]] = 1;
228 std::vector<unsigned int>
tmp;
229 tmp.push_back(Cls[i][c1]);
230 tmp.push_back(Cls[j][c2]);
231 matchedclusters.push_back(tmp);
232 matched[Cls[i][
c1]] = 1;
233 matched[Cls[j][
c2]] = 1;
241 for (
size_t i = 0; i < matchedclusters.size(); ++i) {
242 if (matchedclusters[i].
size())
244 for (
size_t j = 0; j < matchedclusters[i].size(); ++j) {
249 return matchedclusters;
code to link reconstructed objects back to the MC truth information
MaybeLogger_< ELseverityLevel::ELsev_info, true > LogVerbatim
double GetXTicksOffset(int p, int t, int c) const
Float_t x1[n_points_granero]
Planes which measure Z direction.
decltype(auto) constexpr end(T &&obj)
ADL-aware version of std::end.
decltype(auto) constexpr size(T &&obj)
ADL-aware version of std::size.
bool SortByWire(art::Ptr< recob::Hit > const &h1, art::Ptr< recob::Hit > const &h2)
unsigned int NumberTimeSamples() const
double ConvertTicksToX(double ticks, int p, int t, int c) const
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
unsigned int Nplanes(TPCID const &tpcid=tpc_zero) const
Returns the total number of planes in the specified TPC.
decltype(auto) constexpr begin(T &&obj)
ADL-aware version of std::begin.
recob::tracking::Plane Plane
Length_t WirePitch(PlaneID const &planeid=plane_zero) const
Returns the distance between two consecutive wires.