LDMX Software
TrigScintTrackProducer.cxx
1#include "TrigScint/TrigScintTrackProducer.h"
2
3#include <iterator> // std::next
4#include <map>
5
6namespace trigscint {
7
9 maxDelta_ = ps.getParameter<double>(
10 "delta_max"); // max distance to consider adding in a cluster to track
11 seeding_collection_ = ps.getParameter<std::string>(
12 "seeding_collection"); // probably tagger pad, "TriggerPadTagClusters"
13 input_collections_ = ps.getParameter<std::vector<std::string>>(
14 "further_input_collections"); // {"TriggerPadUpClusters" ,
15 // "TriggerPadDnClusters" }
16 output_collection_ = ps.getParameter<std::string>("output_collection");
17 passName_ = ps.getParameter<std::string>("input_pass_name");
18 verbose_ = ps.getParameter<int>("verbosity");
19 vertBarStartIdx_ = ps.getParameter<int>("vertical_bar_start_index");
20 nBarsY_ = ps.getParameter<int>("number_horizontal_bars");
21 barWidth_y_ = ps.getParameter<double>("horizontal_bar_width");
22 barGap_y_ = ps.getParameter<double>("horizontal_bar_gap");
23 nBarsX_ = ps.getParameter<int>("number_vertical_bars");
24 barWidth_x_ = ps.getParameter<double>("vertical_bar_width");
25 barGap_x_ = ps.getParameter<double>("vertical_bar_gap");
26 skipLast_ = ps.getParameter<bool>("allow_skip_last_collection");
27
28 // TO DO: allow any number of input collections
29
30 if (verbose_) {
31 ldmx_log(info) << "In TrigScintTrackProducer: configure done!" << std::endl;
32 ldmx_log(info) << "Got parameters: \nSeeding: " << seeding_collection_
33 << "\nTolerance: " << maxDelta_
34 << "\nInput: " << input_collections_.at(0) << " and "
35 << input_collections_.at(1)
36 << "\nInput pass name: " << passName_
37 << "\nAllow tracks with no hit in last collection: "
38 << skipLast_
39 << "\nVertical bar start index: " << vertBarStartIdx_
40 << "\nNumber of horizontal bars: " << nBarsY_
41 << "\nHorizontal bar width: " << barWidth_y_
42 << "\nHorizontal bar gap: " << barGap_y_
43 << "\nNumber of vertical bars: " << nBarsX_
44 << "\nVertical bar width: " << barWidth_x_
45 << "\nVertical bar gap: " << barGap_x_
46 << "\nOutput: " << output_collection_
47 << "\nVerbosity: " << verbose_;
48 }
49
50 yConvFactor_ =
51 (barWidth_y_ + barGap_y_) /
52 2.; // each bar only goes half this distance up (overlap/zig-zag)
53 yStart_ = -(nBarsY_ * (barWidth_y_ + barGap_y_) - barGap_y_) /
54 2.; // half height of pad
55 xConvFactor_ =
56 barWidth_x_ +
57 barGap_x_; // each bar goes entire distance sideways (no overlap)
58 xStart_ = -(nBarsX_ * (barWidth_x_ + barGap_x_) - barGap_x_) /
59 2.; // half width of pad
60
61 return;
62}
63
65 // parameters.
66 // one pad cluster collection to use as seed
67 // a vector with the other two
68 // a maximum distance between seed centroid and other pad clusters
69 // allowed to belong to the same track
70 // an output collection name a verbosity controller
71
72 if (verbose_) {
73 ldmx_log(debug)
74 << "TrigScintTrackProducer: produce() starts! Event number: "
75 << event.getEventHeader().getEventNumber();
76 }
77 if (!event.exists(seeding_collection_, passName_)) {
78 ldmx_log(info) << "No collection called " << seeding_collection_
79 << "; skipping event";
80 // << "; still, not skipping event";
81 return;
82 }
83
84 if (!event.exists(seeding_collection_, passName_)) {
85 ldmx_log(info) << "No collection called " << seeding_collection_
86 << "; skipping event";
87 return;
88 }
89 const auto seeds{event.getCollection<ldmx::TrigScintCluster>(
90 seeding_collection_, passName_)};
91 uint numSeeds = seeds.size();
92
93 if (verbose_) {
94 ldmx_log(debug) << "Got track seeding cluster collection "
95 << seeding_collection_ << " with " << numSeeds
96 << " entries ";
97 }
98
99 if (!event.exists(input_collections_.at(0), passName_)) {
100 ldmx_log(info) << "No collection called " << input_collections_.at(0)
101 << "; skipping event";
102 // << "; still, not skipping event";
103
104 return;
105 }
106 const auto clusters_pad1{event.getCollection<ldmx::TrigScintCluster>(
107 input_collections_.at(0), passName_)};
108
109 if (!event.exists(input_collections_.at(1), passName_)) {
110 ldmx_log(info) << "No collection called "
111 << input_collections_.at(1)
112 // << "; still, not skipping event";
113 << "; skipping event";
114 std::vector<ldmx::TrigScintTrack> empty{};
115 event.add(output_collection_, empty);
116 return;
117 }
118
119 const auto clusters_pad2{event.getCollection<ldmx::TrigScintCluster>(
120 input_collections_.at(1), passName_)};
121
122 if (verbose_) {
123 ldmx_log(debug) << "Got the other two pad collections:"
124 << input_collections_.at(0) << " with "
125 << clusters_pad1.size() << " entries, and "
126 << input_collections_.at(1) << " with "
127 << clusters_pad2.size() << " entries.";
128 }
129
130 std::vector<ldmx::TrigScintTrack> cleanedTracks;
131 std::vector<ldmx::TrigScintTrack> cleanedTracksY;
132 std::vector<ldmx::TrigScintTrack> cleanedTracksX;
133
134 // loop over the clusters in the seeding pad collection, if there are clusters
135 // in all pads
136 // bool skipDn = false;
137 if (numSeeds && clusters_pad1.size()) {
138 // could check this explicitly here: and then just get out of all checks on
139 // the dn pad immediately
140 // if (! clusters_pad2.size())
141 // skipDn = true ;
142 for (const auto &seed : seeds) {
143 // for each seed, search through the other two pads to match all clusters
144 // with centroids within tolerance to tracks
145 float centroid = seed.getCentroid();
146
147 std::vector<ldmx::TrigScintTrack> trackCandidates;
148
149 if (verbose_ > 1) {
150 ldmx_log(debug) << "Got seed with centroid " << centroid;
151 }
152
153 // reset for each seed
154 // bool madeTrack = false;
155
156 for (const auto &cluster1 : clusters_pad1) {
157 if (verbose_ > 1) {
158 ldmx_log(debug) << "\tGot pad1 cluster with centroid "
159 << cluster1.getCentroid();
160 }
161 if (fabs(cluster1.getCentroid() - centroid) < maxDelta_ ||
162 (centroid >= vertBarStartIdx_ // then in vertical bars
163 && seed.getCentroidX() == cluster1.getCentroidX())) {
164 // use geometry y overlap scheme to see if this is really a match in x
165 // should be done in a map
166
167 if (centroid >= vertBarStartIdx_ &&
168 seed.getCentroidY() < cluster1.getCentroidY()) {
169 // impossible combination
170 if (verbose_ > 1) {
171 ldmx_log(debug) << "\tSkipping impossible x cluster combination "
172 "with y flags (tag up) ("
173 << seed.getCentroidY() << " "
174 << cluster1.getCentroidY() << ")";
175 }
176 continue;
177 }
178
179 // else: first (possible) match! loop through next pad too
180
181 if (verbose_ > 1) {
182 ldmx_log(debug) << "\t\tIt is close enough!. Check pad2";
183 }
184
185 // try making third pad clusters an optional part of track
186
187 std::vector<ldmx::TrigScintCluster> clusterVec = {seed, cluster1};
188
189 bool hasMatchDn = false;
190
191 for (const auto &cluster2 : clusters_pad2) {
192 if (verbose_ > 1) {
193 ldmx_log(debug) << "\tGot pad2 cluster with centroid "
194 << cluster2.getCentroid();
195 }
196
197 if (fabs(cluster2.getCentroid() - centroid) < maxDelta_ ||
198 (centroid >= vertBarStartIdx_ // then in vertical bars
199 && seed.getCentroidX() == cluster2.getCentroidX())) {
200 // use geometry y overlap scheme to see if this is really a match
201 // in x
202
203 if (centroid >= vertBarStartIdx_ &&
204 (seed.getCentroidY() < cluster2.getCentroidY() ||
205 cluster1.getCentroidY() >
206 cluster2.getCentroidY())) { // impossible
207 if (verbose_ > 1) {
208 ldmx_log(debug)
209 << "\tSkipping impossible x cluster combination with y "
210 "flags (tag up dn) ("
211 << seed.getCentroidY() << " " << cluster1.getCentroidY()
212 << " " << cluster2.getCentroidY() << ")";
213 }
214 continue;
215 }
216
217 // first match! loop through next pad too
218
219 if (verbose_ > 1) {
220 ldmx_log(debug) << "\t\tIt is close enough!. Make a track";
221 }
222
223 // only make this vector now! this ensures against hanging
224 // clusters with indices from earlier in the loop
225 std::vector<ldmx::TrigScintCluster> threeClusterVec = {
226 seed, cluster1, cluster2};
227
228 /*
229 // here we could break if we didn't want to allow all possible
230 combinations madeTrack=true; break; //we're done with this
231 iteration once there's a track made
232 */
233 // make a track
234 ldmx::TrigScintTrack track = makeTrack(threeClusterVec);
235 trackCandidates.push_back(track);
236 hasMatchDn = true;
237 } // if match in pad2
238 } // over clusters in pad2
239 // if there was no match to this in pad 2, make a track with just
240 // these two clusters
241 if (!hasMatchDn &&
242 skipLast_) { // we allow skipping last pad if needed
243 ldmx::TrigScintTrack track = makeTrack(clusterVec);
244 trackCandidates.push_back(track);
245 }
246
247 } // if possible (x,)y match in pad1
248 /*
249 //same here
250 if (madeTrack)
251 break;
252 */
253
254 } // over clusters in pad1
255
256 // continue to next seed if 0 track candidates
257 if (trackCandidates.size() == 0) continue;
258
259 int keepIdx = 0;
260 float minResidual = 1000; // some large number
261
262 // no need to choose between only one candidate track
263 if (trackCandidates.size() > 1) {
264 // now for each seed, pick only the track with the smallest residual.
265
266 if (verbose_) {
267 ldmx_log(debug) << "Got " << trackCandidates.size()
268 << " tracks to check.";
269 }
270
271 for (uint idx = 0; idx < trackCandidates.size(); idx++) {
272 if ((trackCandidates.at(idx)).getResidual() < minResidual) {
273 keepIdx = (int)idx;
274 minResidual =
275 (trackCandidates.at(idx)).getResidual(); // update minimum
276
277 if (verbose_ > 1) {
278 ldmx_log(debug)
279 << "Track at index " << idx
280 << " has smallest residual so far: " << minResidual;
281 }
282
283 } // finding min residual
284 } // over track candidates
285 } // if more than 1 to choose from
286
287 // store the track at keepIdx, if there was one we made it this far and
288 // keepIdx is 0 or has been updated to the smallest residual track idx
289 // if (keepIdx >= 0) {
290 tracks_.push_back(trackCandidates.at(keepIdx));
291 if (verbose_) {
292 ldmx_log(debug) << "Kept track at index " << keepIdx;
293 (trackCandidates.at(keepIdx)).Print();
294 }
295 //}
296 } // over seeds
297
298 // done here if there were no tracks found
299 if (tracks_.size() == 0) {
300 if (verbose_) {
301 ldmx_log(debug) << "No tracks found!";
302 }
303 std::vector<ldmx::TrigScintTrack> empty{};
304 event.add(output_collection_, empty);
305 return;
306 }
307 // now, if there are multiple seeds sharing the same downstream hits, this
308 // should also be remedied with a selection on min residual.
309
310 // The logic of this loop kind of assumes I can remove tracks immediately --
311 // that way I can do pairwise checks between more tracks within a single
312 // loop. But for now I haven't figured out how to erase elements in a fool
313 // proof way. So I iterate over a vector...
314
315 std::vector keepIndices(tracks_.size(), 1);
316 if (verbose_ > 1)
317 ldmx_log(debug) << "vector of indices to keep has size "
318 << keepIndices.size();
319
320 for (uint idx = tracks_.size() - 1; idx > 0; idx--) {
321 // since we start in one end, we only have to check matches in one
322 // direction
323 ldmx::TrigScintTrack track = tracks_.at(idx);
324 for (int idxComp = idx - 1; idxComp >= 0; idxComp--) {
325 if (verbose_ > 1)
326 ldmx_log(debug) << "In track disambiguation loop, idx points at "
327 << idx << " and prev idx points at " << idxComp;
328
329 ldmx::TrigScintTrack nextTrack = tracks_.at(idxComp);
330
331 // no need to start pulling constituents from tracks that are
332 // ridiculously far apart
333 if (fabs(track.getCentroid() - nextTrack.getCentroid() <
334 3 * maxDelta_)) {
335 std::vector<ldmx::TrigScintCluster> consts_1 =
336 track.getConstituents();
337 std::vector<ldmx::TrigScintCluster> consts_2 =
338 nextTrack.getConstituents();
339 if (verbose_ > 1)
340 ldmx_log(debug)
341 << "In track disambiguation loop, got the two tracks, "
342 "with nConstituents "
343 << consts_1.size() << " and " << consts_2.size()
344 << ", respectively. ";
345 // let's do "if either cluster is shared" right now... but could also
346 // have it settable to use a stricter cut: an AND
347 if (consts_1[1].getCentroid() == consts_2[1].getCentroid() ||
348 (consts_1.size() > 2 && consts_2.size() > 2 &&
349 consts_1[2].getCentroid() ==
350 consts_2[2]
351 .getCentroid())) { // we have overlap downstream of the
352 // seeding pad. probably, one cluster
353 // in seeding pad is noise
354
355 if (verbose_ > 1) {
356 ldmx_log(debug) << "Found overlap! Tracks at index " << idx
357 << " and " << idxComp;
358 (tracks_.at(idx)).Print();
359 (tracks_.at(idxComp)).Print();
360 }
361
362 if ((tracks_.at(idx)).getResidual() <
363 (tracks_.at(idxComp)).getResidual()) {
364 // next track (lower index) is a worse choice, remove its flag for
365 // keeping
366 keepIndices.at(idxComp) = 0;
367 } else // prefer next track over current. remove current track's
368 // keep
369 // flag
370 keepIndices.at(idx) = 0;
371 /*}
372 else {
373 tracks_.erase(itNext);
374 // removeIdx.push_back(idx+1);
375 // we might see the same index two times in the loop in this
376 case, if there are three seeds sharing the same clusters
377 downstream.
378 // then the third only gets removed if it's even worse than
379 the second.
380 // one could deal with this with an extra overlap check. not
381 sure we will be in this situation any time soon though.
382 }*/
383 } // over matching/overlapping tracks
384 } // over tracks close enough to share constituents
385 } // over constructed tracks at other indices, to match
386 } // over constructed tracks
387
388 for (uint idx = 0; idx < tracks_.size(); idx++) {
389 if (verbose_ > 1) {
390 ldmx_log(debug) << "keep flag for idx " << idx << " is "
391 << keepIndices.at(idx);
392 }
393 if (keepIndices.at(idx)) { // this hasn't been flagged for removal
394
395 cleanedTracks.push_back(tracks_.at(idx));
396
397 if (verbose_) {
398 ldmx_log(debug) << "After cleaning, keeping track at index " << idx
399 << ": Centroid = " << (tracks_.at(idx)).getCentroid()
400 << "; CentroidX = "
401 << (tracks_.at(idx)).getCentroidX()
402 << "; CentroidY = "
403 << (tracks_.at(idx)).getCentroidY()
404 << "; track PE = " << (tracks_.at(idx)).getPE();
405 // (tracks_.at(idx)).Print();
406 }
407 } // if index flagged for keeping
408 } // over all (uniquely seeded) tracks in the event
409 /*
410 if (verbose_ ) {
411 for (uint idx=0; idx < tracks_.size(); idx++){
412 ldmx_log(debug)<< "Keeping track at index " << idx << ":";
413 (tracks_.at(idx)).Print();
414 }
415 }
416 */
417
418 if (verbose_) {
419 ldmx_log(debug) << "Running track x,y matching ";
420 }
421
422 if (cleanedTracks.size() > 0) {
423 matchXYTracks(cleanedTracks);
424 std::vector<ldmx::TrigScintTrack> matchedTracks =
425 cleanedTracks; // don't know why this copying needs to happen but it
426 // does
427 // std::vector<ldmx::TrigScintTrack> matchXYTracks( cleanedTracks
428 //); std::vector<ldmx::TrigScintTrack> matchedTracks =
429 // matchXYTracks( cleanedTracks );
430 for (auto trk : matchedTracks) {
431 /* for (uint idx = 0; idx < tracks_.size(); idx++) {
432 if (verbose_ > 1) {
433 ldmx_log(debug) << "keep flag for idx " << idx << " is "
434 << keepIndices.at(idx);
435 }
436 if (keepIndices.at(idx)) { // this hasn't
437 been flagged for removal
438 //check if channel nb is above that of horizontal bars
439 if (tracks_.at(idx).getCentroid() >= vertBarStartIdx_)
440 */
441 if (trk.getCentroid() >= vertBarStartIdx_)
442 cleanedTracksX.push_back(trk); // acks_.at(idx));
443 else
444 cleanedTracksY.push_back(trk); // acks_.at(idx));
445 // cleanedTracksY.push_back(trk);
446 if (verbose_ > 1) {
447 float centr = trk.getCentroid(); // tracks_.at(idx).getCentroid(); //
448 std::string collStr = centr >= vertBarStartIdx_ ? "X" : "Y";
449 collStr = output_collection_ + collStr;
450 ldmx_log(debug) << "saving track with centroid " << centr
451 << " to output track collection " << collStr;
452 }
453 // }
454 }
455 }
456
457 } // if there are clusters in all pads
458 else if (verbose_) {
459 ldmx_log(info)
460 << "Not all pads had clusters; (maybe) skipping tracking attempt";
461 }
462
463 if (verbose_) {
464 ldmx_log(debug) << "Done with tracking step. ";
465 }
466
467 event.add(output_collection_, cleanedTracks);
468 // event.add(output_collection_, matchedTracks);
469
470 event.add(output_collection_ + "Y", cleanedTracksY);
471 event.add(output_collection_ + "X", cleanedTracksX);
472
473 tracks_.resize(0);
474
475 return;
476}
477
478ldmx::TrigScintTrack TrigScintTrackProducer::makeTrack(
479 std::vector<ldmx::TrigScintCluster> clusters) {
480 // for now let's keep a straight, unweighted centroid
481 // consider the possibility that at least one cluster has a centroid
482 // identically == 0. then we need to shift them by 1 if we want to do energy
483 // weighted track centroid later. but no need now
485 float centroid = 0;
486 float centroidX = 0;
487 float centroidY = 0;
488 float beamEfrac = 0;
489 float pe = 0;
490 for (uint i = 0; i < clusters.size(); i++) {
491 centroid += (clusters.at(i)).getCentroid();
492 centroidX += (clusters.at(i)).getCentroidX();
493 centroidY += (clusters.at(i)).getCentroidY();
494 tr.addConstituent(clusters.at(i));
495 beamEfrac += (clusters.at(i)).getBeamEfrac();
496 pe += (clusters.at(i)).getPE();
497 }
498 centroid /= clusters.size();
499 centroidX /= clusters.size();
500 if (centroid >= vertBarStartIdx_) {
501 if (verbose_) {
502 ldmx_log(debug)
503 << " -- In makeTrack made vertical bar track with centroid "
504 << centroid << " and y flag sum " << centroidY;
505 // try commenting this to check if that helps with an out-of-bounds
506 // problem
507 // << " from clusters with y centroids";
508 // for (uint i = 0; i < clusters.size(); i++)
509 // ldmx_log(debug) << "\tpad " << i << ": centroidY "
510 // << (clusters.at(i)).getCentroidY();
511 }
512 // then the sum of centroid y is 0, 2, 4 or 6
513 // we have 4 divisions, so, the center of it should be divNb/8
514 // (or rather, that's where channel nBars/8 begins)
515 // and then a factor 2 for the zig-zag pattern
516 centroidY = (centroidY + 1) * 2 * nBarsY_ / 8.;
517 // TODO: here we could instead just use quadrant indices 0-3 by dividing by
518 // 2 but that would mean that in the raw, x and y track centroidY would mean
519 // different things
520 if (verbose_) ldmx_log(debug) << " -- new centroidY = " << centroidY;
521 } else
522 centroidY /= clusters.size();
523
524 beamEfrac /= clusters.size();
525 pe /= clusters.size();
526
527 float residual = 0;
528 for (uint i = 0; i < clusters.size(); i++)
529 residual += ((clusters.at(i)).getCentroid() - centroid) *
530 ((clusters.at(i)).getCentroid() - centroid);
531 residual = sqrt(residual / clusters.size());
532
533 tr.setCentroid(centroid);
534 tr.setCentroidX(centroidX);
535 tr.setCentroidY(centroidY);
536 tr.setResidual(residual);
537 tr.setBeamEfrac(beamEfrac);
538 tr.setPE(pe);
539
540 if (verbose_) {
541 ldmx_log(debug) << " -- In makeTrack made track with centroid "
542 << centroid << " and residual " << residual << " and pe "
543 << pe << " from clusters with centroids";
544 for (uint i = 0; i < clusters.size(); i++)
545 ldmx_log(debug) << "\tpad " << i << ": centroid "
546 << (clusters.at(i)).getCentroid();
547 }
548
549 return tr;
550}
551
552// std::vector<ldmx::TrigScintTrack> TrigScintTrackProducer::matchXYTracks(
553void TrigScintTrackProducer::matchXYTracks(
554 std::vector<ldmx::TrigScintTrack> &tracks) {
555 // map quadrant nb to track (can be multiple per quadrant)
556 std::multimap<int, int>
557 yIdxQuadMap; // key = quad, val = track index in collection
558 std::multimap<int, int> xIdxQuadMap;
559
560 std::multimap<int, ldmx::TrigScintTrack> yQuadMap;
561 std::multimap<int, ldmx::TrigScintTrack> xQuadMap;
562 // map track in quadrant back to index in entire track collection
563 // used for updating collection track variables
564 std::map<ldmx::TrigScintTrack, int> yTrackMap;
565 std::map<ldmx::TrigScintTrack, int> xTrackMap;
566
567 uint trkIdx = -1;
568 for (auto trk : tracks) {
569 trkIdx++;
570 // 1. get the y bar tracks with centroidX = -1
571 if (trk.getCentroidX() == -1) {
572 if (verbose_)
573 ldmx_log(debug) << " -- In matchXYTracks found y track at "
574 << trk.getCentroidY() << "; mapping to quad "
575 << (int)trk.getCentroidY() / 8 << " with trk index "
576 << trkIdx;
577 // 2. order them... or map them to quadrants. note that there are 2 layers
578 // so 2*nBarsY_/4 channels per quadrant
579 yQuadMap.insert(std::make_pair((int)(trk.getCentroidY() / 8), trk));
580 yTrackMap[trk] = trkIdx;
581 yIdxQuadMap.insert(std::make_pair((int)(trk.getCentroidY() / 8), trkIdx));
582
583 } else { // 3. get the remaining tracks (from vertical bars) and map them
584 // (back) to (middle of) quadrants
585 xQuadMap.insert(std::make_pair((int)(trk.getCentroidY() / 8), trk));
586 xTrackMap[trk] = trkIdx;
587 xIdxQuadMap.insert(std::make_pair((int)(trk.getCentroidY() / 8), trkIdx));
588 if (verbose_)
589 ldmx_log(debug) << " -- In matchXYTracks found x track at (x,y) = ("
590 << trk.getCentroidX() << ", " << trk.getCentroidY()
591 << "); mapping to quad " << (int)trk.getCentroidY() / 8
592 << " with trk index " << trkIdx;
593 }
594 }
595
596 // 4a
597 //
598 // 1) here use the geometry? if we can assume perfect alignment we can take
599 // width and nBars and take nBars/2 as origin
600 // --- now do the matching ---
601
602 // if there is no useful matching to be done: these are the pad width wide
603 // numbers
604 float x0 = 0;
605 float sx0 = fabs(xStart_); // this should be half the pad... could also set
606 // it to full beam spot width
607 float sy0 = fabs(yStart_) /
608 4.; // yStart_ is half the pad, so this should be half a quadrant
609
610 // assume at least one y track. will have to figure out if there is ever a
611 // reason to use an isolated x track in its place.
612 for (auto yitr = yQuadMap.begin(); yitr != yQuadMap.end(); ++yitr) {
613 int nYinQuad = yQuadMap.count((*yitr).first);
614 int nXinQuad = xQuadMap.count((*yitr).first);
615 float y{-9999.}, sy{-9999.}, x{-9999.}, x1{-9999.}, x2{-9999.}, sx1{-9999.},
616 sx2{-9999.}, y1{-9999.}, y2{-9999.}, sy1{-9999.}, sy2{-9999.};
617 // quad midpoint:
618 float y0 = (*yitr).first * 8 + sy0;
619 float sx =
620 1. / 2 * xConvFactor_; // rely on x precision being one single bar
621 // width; always used unless x is undeterminable
622
623 // check all x first
624 // do the easiest first:
625 if (nXinQuad == 0) { // then there's no hope of setting a better x here
626 // just use the beam spot width... and center of pad
627 x = x0;
628 sx = sx0;
629 if (verbose_)
630 ldmx_log(debug) << "\t\t\t no x info in quad " << (*yitr).first
631 << "; will set x to middle of pad, pad half-width as "
632 "precision: set (x, sx)=("
633 << x << ", " << sx << ")";
634 } // 0 x tracks in quadrant
635 else if (nXinQuad ==
636 1) { // slightly harder: 1 x track -- might be easy if
637 // it's just one y track; if several, need to
638 // think about overlaps. but in overlap case, just
639 // revert to setting x0 and sx0, when we know
640 auto xitr = xQuadMap.find((*yitr).first);
641 x = ((*xitr).second).getCentroidX() * xConvFactor_ + xStart_;
642
643 if (verbose_)
644 ldmx_log(debug) << "\t\t\t 1 x in quad " << (*yitr).first
645 << ", getting (x, sx)=(" << x << ", " << sx << ")";
646 } // 1 x track in quadrant
647 else if (nXinQuad == 2) { // finally if we have two tracks, get x1 and x2
648 // and decide later how to use them
649 // don't think we want to experiment with discerning three overlapping
650 // tracks, so not >= 2
651 // continue; //debugging: skip for now -- didn't help
652 auto xitr1 = xQuadMap.lower_bound((*yitr).first);
653 auto xitr2 = xQuadMap.upper_bound((*yitr).first);
654 xitr2--; // upper_bound points to next element
655
656 if (xitr1 != xitr2) { // should be true already but...
657 x1 = ((*xitr1).second).getCentroidX() * xConvFactor_ + xStart_;
658 x2 = ((*xitr2).second).getCentroidX() * xConvFactor_ + xStart_;
659 sx1 = xConvFactor_ / 2.; // 1 bar width
660 sx2 = sx1;
661 x = (x1 + x2) / 2.;
662 sx = fabs(x1 - x2) / 2 *
663 xConvFactor_; // rely on x precision being one single pad width
664 if (verbose_)
665 ldmx_log(debug) << "\t\t -- 2 x in quad: setting y track x "
666 "coordinate to midpoint";
667 }
668 } // if 2 x tracks in quad
669
670 // ok! over y:
671 // can skip 0 y case by construction
672 if (nYinQuad == 1) { // we can already now tell what the y coordinate and
673 // its precision is
674 y = ((*yitr).second).getCentroidY() * yConvFactor_ + yStart_;
675 sy = ((*yitr).second).getResidual() * yConvFactor_;
676 if (sy == 0)
677 sy = 1. / 2 * yConvFactor_; // if all clusters lined up, assign
678 // precision of 1 bar width
679
680 if (nXinQuad <= 1) {
681 // 4. every quadrant which just has one of each --> done ;
682 // b) set the sx, sy of the x track now, using the residuals from the
683 // other b1) special case: no x tracks; then x, sx have been set above
684 auto xidx = xIdxQuadMap.find((*yitr).first);
685 auto yidx = yIdxQuadMap.find((*yitr).first);
686 tracks.at((*xidx).second).setPosition(x, y);
687 tracks.at((*xidx).second).setSigmaXY(sx, sy);
688 tracks.at((*yidx).second).setPosition(x, y);
689 tracks.at((*yidx).second).setSigmaXY(sx, sy);
690 if (verbose_)
691 ldmx_log(debug) << "\t\t\t in quad " << (*yitr).first
692 << ", set (x, y) = (" << x << ", " << y
693 << ") and (sx, sy) = " << sx << ", " << sy << ")";
694 continue;
695 }
696 } // 1 y, 0 or 1 x track in quadrant
697
698 if (verbose_)
699 ldmx_log(debug) << "\t\t in quad " << (*yitr).first
700 << ", not single x,y tracks: " << nXinQuad << " of x and "
701 << nYinQuad << " of y";
702
703 if (nYinQuad == 2) { // let's start here and see if we can do >= 2 later
704 // here one could do sth to avoid checking the other y track again in the
705 // outermost loop over y
706 auto yitr1 = yQuadMap.lower_bound((*yitr).first);
707 auto yitr2 = yQuadMap.upper_bound((*yitr).first);
708 yitr2--; // back up once
709 y1 = ((*yitr1).second).getCentroidY() * yConvFactor_ + yStart_;
710 y2 = ((*yitr2).second).getCentroidY() * yConvFactor_ + yStart_;
711 sy1 = ((*yitr1).second).getResidual() * yConvFactor_;
712 sy2 = ((*yitr2).second).getResidual() * yConvFactor_;
713 y = (y1 + y2) / 2.;
714 sy = fabs(y1 - y2) / 2 * yConvFactor_;
715 if (verbose_)
716 ldmx_log(debug)
717 << "\t\t -- 2 y in quad: setting x track y coordinate to midpoint";
718 } // 2y in quad
719
720 if (nYinQuad == 1 &&
721 nXinQuad == 2) { // don't think we want to experiment with discerning
722 // three overlapping tracks, so not >= 2
723
724 // first: set the y track coordinates to x = the mid of x tracks, y = y
725 // of y track
726 auto yidx = yIdxQuadMap.find((*yitr).first);
727 tracks.at((*yidx).second).setPosition(x, y);
728 tracks.at((*yidx).second).setSigmaXY(sx, sy);
729
730 int minOverlapPE_ = 250;
731 if (((*yitr).second).getPE() < minOverlapPE_) {
732 // can't tell, really, that either of these belong to the y track. so.
733 // let them keep their own x coordinate but set y to quadrant midpoint,
734 // with uncertainty +/- half quadrant width (1/8 of pad height)
735 y = y0;
736 sy = sy0;
737 if (verbose_)
738 ldmx_log(debug) << "\t\t -- Can't tell which x track should be "
739 "matched to single y track. Setting both x track "
740 "coordinates to y quadrant value:";
741 } // if can't assume overlap
742 else if (verbose_)
743 ldmx_log(debug) << "\t\t -- Found large PE count ("
744 << ((*yitr).second).getPE() << " > " << minOverlapPE_
745 << "), suggesting overlap! Setting both x track "
746 "coordinates to y track value:";
747
748 // consider making two x tracks out if this one, and, anyway have to set
749 // their average as the y track x cocordinate
750 // EXPERIMENTAL : apply only to x tracks, which can be disregarded for
751 // electron counting
752 if (verbose_)
753 ldmx_log(debug) << "\t\t -- (x1, x2, y) = (" << x1 << ", " << x2
754 << ", " << y << ") and (sx1, sx2, sy) = " << sx1 << ", "
755 << sx2 << ", " << sy << ")";
756
757 // now set x track coordinates according to overlap check result
758 auto xidx1 = xIdxQuadMap.lower_bound((*yitr).first);
759 auto xidx2 = xIdxQuadMap.upper_bound((*yitr).first);
760 xidx2--; // upper_bound points to (last+1) element
761 tracks.at((*xidx1).second).setPosition(x1, y);
762 tracks.at((*xidx1).second).setSigmaXY(sx1, sy);
763 tracks.at((*xidx2).second).setPosition(x2, y);
764 tracks.at((*xidx2).second).setSigmaXY(sx2, sy);
765
766 } // 1 y, 2 x tracks in the quadrant
767 else if (nYinQuad == 2 && nXinQuad == 1) {
768 // 5b) if there are more y than x: could be an overlap
769
770 // first: set the x track coordinates to x = x of x track, y = the mid of
771 // y tracks
772 auto xidx = xIdxQuadMap.find((*yitr).first);
773 tracks.at((*xidx).second).setPosition(x, y);
774 tracks.at((*xidx).second).setSigmaXY(sx, sy);
775
776 auto xitr = xQuadMap.lower_bound((*yitr).first);
777 int minOverlapPE_ = 300;
778 if (((*xitr).second).getPE() < minOverlapPE_) {
779 if (verbose_)
780 ldmx_log(debug)
781 << "\t\t just 1 x track with not-unusual PE in the quad -- can't "
782 "match; setting mid-point values for x ";
783 x = x0;
784 sx = sx0;
785 } // if can't assume overlap
786 else {
787 // consider making two x tracks out if this one, and, anyway have to set
788 // their average as the y track x cocordinate
789 // EXPERIMENTAL : apply only to x tracks, which can be disregarded for
790 // electron counting
791 if (verbose_)
792 ldmx_log(debug) << "\t\t -- Found large PE count ("
793 << ((*xitr).second).getPE() << " > " << minOverlapPE_
794 << ") in x track, suggesting overlap! Setting both y "
795 "track coordinates to x track value:";
796 } // if can assume overlap
797 if (verbose_)
798 ldmx_log(debug) << "\t\t -- (x, y1, y2) = (" << x << ", " << y1 << ", "
799 << y2 << ") and (sx, sy1, sy2) = " << sx << ", " << sy1
800 << ", " << sy2 << ")";
801
802 auto yidx1 = yIdxQuadMap.lower_bound((*yitr).first);
803 auto yidx2 = yIdxQuadMap.upper_bound((*yitr).first);
804 yidx2--; // upper_bound points to next element
805 tracks.at((*yidx1).second).setPosition(x, y1);
806 tracks.at((*yidx1).second).setSigmaXY(sx, sy1);
807 tracks.at((*yidx2).second).setPosition(x, y2);
808 tracks.at((*yidx2).second).setSigmaXY(sx, sy2);
809
810 } // 2 y and 1 x track in quad
811 else if (nYinQuad == 2 && nXinQuad == 2) {
812 // MIDPONTS ALL OVER!
813 auto xidx1 = xIdxQuadMap.lower_bound((*yitr).first);
814 auto xidx2 = xIdxQuadMap.upper_bound((*yitr).first);
815 xidx2--;
816 auto yidx1 = yIdxQuadMap.lower_bound((*yitr).first);
817 auto yidx2 = yIdxQuadMap.upper_bound((*yitr).first);
818 yidx2--;
819
820 if (yIdxQuadMap.find((*yitr).first) == yIdxQuadMap.end())
821 ldmx_log(error) << "The two y tracks in the same quadrant at "
822 << (*yitr).first
823 << " appear to not be found in the y track map! "
824 "investigate. Note that yidx1.first = "
825 << (*yidx1).first
826 << " and yidx2.first = " << (*yidx2).first;
827 else {
828 tracks.at((*xidx1).second).setPosition(x1, y);
829 tracks.at((*xidx1).second).setSigmaXY(sx1, sy);
830 tracks.at((*xidx2).second).setPosition(x2, y);
831 tracks.at((*xidx2).second).setSigmaXY(sx2, sy);
832
833 tracks.at((*yidx1).second).setPosition(x, y1);
834 tracks.at((*yidx1).second).setSigmaXY(sx, sy1);
835 tracks.at((*yidx2).second).setPosition(x, y2);
836 tracks.at((*yidx2).second).setSigmaXY(sx, sy2);
837
838 if (verbose_)
839 ldmx_log(debug) << "\t\t -- in a 2 x 2 situaiton; midpoint y: " << y
840 << " for both x tracks, midpoint x: " << x
841 << " for both y tracks";
842 }
843 } // if 2 y, 2 x tracks
844
845 if (nXinQuad > 2) {
846 if (verbose_)
847 ldmx_log(debug) << "\t\t -*-*-*- more than 2 x tracks in the same quad "
848 "-- nothing done about the x,y coordinates in this "
849 "situation -- implement if needed!!";
850 }
851 if (nYinQuad > 2) {
852 if (verbose_)
853 ldmx_log(debug) << "\t\t -*-*-*- more than 2 y tracks in the same quad "
854 "-- nothing done about the x,y coordinates in this "
855 "situation -- implement if needed!!";
856 }
857
858 } // over y tracks
859
860 yQuadMap.clear();
861 xQuadMap.clear();
862
863 // return tracks;
864}
865
867 ldmx_log(debug) << "Process starts!";
868
869 return;
870}
871
873 ldmx_log(debug) << "Process ends!";
874
875 return;
876}
877
878} // namespace trigscint
879
880DECLARE_PRODUCER_NS(trigscint, TrigScintTrackProducer);
#define DECLARE_PRODUCER_NS(NS, CLASS)
Macro which allows the framework to construct a producer given its name during configuration.
Implements an event buffer system for storing event data.
Definition Event.h:41
bool exists(const std::string &name, const std::string &passName="", bool unique=true) const
Check for the existence of an object or collection with the given name and pass name in the event.
Definition Event.cxx:92
Class encapsulating parameters for configuring a processor.
Definition Parameters.h:27
T getParameter(const std::string &name) const
Retrieve the parameter of the given name.
Definition Parameters.h:89
Stores cluster information from the trigger scintillator pads.
Represents a track of trigger scintillator clusters.
void setCentroidX(float centroid)
Set the x centroid of the track.
void setResidual(float resid)
Set the detector ID residual of the track.
void setCentroidY(float centroid)
Set the y centroid of the track.
void setPE(float pe)
Set the average cluster pe of the track.
float getCentroid() const
Get the detector ID centroid of the track.
void addConstituent(TrigScintCluster cl)
Add a cluster to the list of track constituents.
void setCentroid(float centroid)
Set the detector ID centroid of the track.
void setBeamEfrac(float e)
Set beam energy fraction of hit.
std::vector< ldmx::TrigScintCluster > getConstituents() const
Get the cluster constituents of the track.
void configure(framework::config::Parameters &ps) override
Callback for the EventProcessor to configure itself from the given set of parameters.
void onProcessEnd() override
Callback for the EventProcessor to take any necessary action when the processing of events finishes,...
void produce(framework::Event &event) override
Process the event and put new data products into it.
void onProcessStart() override
Callback for the EventProcessor to take any necessary action when the processing of events starts,...