LArSoft  v10_04_05
Liquid Argon Software toolkit - https://larsoft.org/
BeamFlashTrackMatchTaggerAlg.cxx
Go to the documentation of this file.
1 
13 
14 #include "fhiclcpp/ParameterSet.h"
15 
23 
24 #include <limits>
25 #include <utility>
26 #include <vector>
27 
28 #include "TH1F.h"
29 #include "TTree.h"
30 #include "TVector3.h"
31 
33  : COSMIC_TYPE_FLASHMATCH(anab::CosmicTagID_t::kFlash_BeamIncompatible)
34  , COSMIC_TYPE_OUTSIDEDRIFT(anab::CosmicTagID_t::kOutsideDrift_Partial)
35  , DEBUG_FLAG(p.get<bool>("RunDebugMode", false))
36  , fMinTrackLength(p.get<float>("MinTrackLength"))
37  , fMinOpHitPE(p.get<float>("MinOpHitPE", 0.1))
38  , fMIPdQdx(p.get<float>("MIPdQdx", 2.1))
39  , fOpDetSaturation(p.get<float>("OpDetSaturation", 200.))
40  , fSingleChannelCut(p.get<float>("SingleChannelCut"))
41  , fCumulativeChannelThreshold(p.get<float>("CumulativeChannelThreshold"))
42  , fCumulativeChannelCut(p.get<unsigned int>("CumulativeChannelCut"))
43  , fIntegralCut(p.get<float>("IntegralCut"))
44  , fMakeOutsideDriftTags(p.get<bool>("MakeOutsideDriftTags", false))
45  , fNormalizeHypothesisToFlash(p.get<bool>("NormalizeHypothesisToFlash"))
46 {}
47 
49  TH1F* hist_flash,
50  TH1F* hist_hyp)
51 {
52  cTree = tree;
53 
54  cOpDetHist_flash = hist_flash;
55  cOpDetHist_flash->SetNameTitle("opdet_hist_flash", "Optical Detector Occupancy, Flash");
56 
57  cOpDetHist_hyp = hist_hyp;
58  cOpDetHist_hyp->SetNameTitle("opdet_hist_hyp", "Optical Detector Occupancy, Hypothesis");
59 
61  cTree->Branch("opdet_hyp", &cOpDetVector_hyp);
62  cTree->Branch("opdet_flash", &cOpDetVector_flash);
63  cTree->Branch("opdet_hist_flash", "TH1F", cOpDetHist_flash);
64  cTree->Branch("opdet_hist_hyp", "TH1F", cOpDetHist_hyp);
65 }
66 
68  std::vector<recob::OpFlash> const& flashVector,
69  std::vector<recob::Track> const& trackVector,
70  std::vector<anab::CosmicTag>& cosmicTagVector,
71  std::vector<size_t>& assnTrackTagVector,
72  Providers_t providers,
74  opdet::OpDigiProperties const& opdigip)
75 {
76  auto const& geom = *providers.get<geo::GeometryCore>();
77  auto const& wireReadoutGeom = art::ServiceHandle<geo::WireReadout const>()->Get();
78 
79  std::vector<const recob::OpFlash*> flashesOnBeamTime;
80  for (auto const& flash : flashVector) {
81  if (!flash.OnBeamTime()) continue;
82  flashesOnBeamTime.push_back(&flash);
83  }
84 
85  //make sure this association vector is initialized properly
86  assnTrackTagVector.resize(trackVector.size(), std::numeric_limits<size_t>::max());
87  cosmicTagVector.reserve(trackVector.size());
88 
89  for (size_t track_i = 0; track_i < trackVector.size(); track_i++) {
90 
91  recob::Track const& track(trackVector[track_i]);
92 
93  if (track.Length() < fMinTrackLength) continue;
94 
95  //get the begin and end points of this track
96  TVector3 const& pt_begin = track.LocationAtPoint<TVector3>(0);
97  TVector3 const& pt_end = track.LocationAtPoint<TVector3>(track.NumberTrajectoryPoints() - 1);
98  std::vector<float> xyz_begin = {(float)pt_begin.x(), (float)pt_begin.y(), (float)pt_begin.z()};
99  std::vector<float> xyz_end = {(float)pt_end.x(), (float)pt_end.y(), (float)pt_end.z()};
100 
101  //check if this track is outside the drift window, and if it is continue
102  if (!InDriftWindow(pt_begin.x(), pt_end.x(), geom)) {
103  if (fMakeOutsideDriftTags) {
104  cosmicTagVector.emplace_back(xyz_begin, xyz_end, 1., COSMIC_TYPE_OUTSIDEDRIFT);
105  assnTrackTagVector[track_i] = cosmicTagVector.size() - 1;
106  }
107  continue;
108  }
109 
110  //get light hypothesis for track
111  std::vector<float> lightHypothesis = GetMIPHypotheses(track, providers, pvs, opdigip);
112 
113  //check compatibility with beam flash
114  bool compatible = false;
115  for (const recob::OpFlash* flashPointer : flashesOnBeamTime) {
116  CompatibilityResultType result =
117  CheckCompatibility(lightHypothesis, flashPointer, geom, wireReadoutGeom);
118  if (result == CompatibilityResultType::kCompatible) compatible = true;
119  if (DEBUG_FLAG) {
120  PrintTrackProperties(track);
121  PrintFlashProperties(*flashPointer);
123  lightHypothesis, flashPointer, geom, wireReadoutGeom, result);
124  }
125  }
126 
127  //make tag
128  float cosmicScore = 1.;
129  if (compatible) cosmicScore = 0.;
130  cosmicTagVector.emplace_back(xyz_begin, xyz_end, cosmicScore, COSMIC_TYPE_FLASHMATCH);
131  assnTrackTagVector[track_i] = cosmicTagVector.size() - 1;
132  }
133 }
134 
135 //this compares the hypothesis to the flash itself.
137  unsigned int const run,
138  unsigned int const event,
139  std::vector<recob::OpFlash> const& flashVector,
140  std::vector<recob::Track> const& trackVector,
141  Providers_t providers,
143  opdet::OpDigiProperties const& opdigip)
144 {
145  auto const& geom = *(providers.get<geo::GeometryCore>());
146  auto const& wireReadoutGeom = art::ServiceHandle<geo::WireReadout const>()->Get();
147 
148  cFlashComparison_p.run = run;
149  cFlashComparison_p.event = event;
150 
151  std::vector<std::pair<unsigned int, const recob::OpFlash*>> flashesOnBeamTime;
152  for (unsigned int i = 0; i < flashVector.size(); i++) {
153  recob::OpFlash const& flash = flashVector[i];
154  if (!flash.OnBeamTime()) continue;
155  flashesOnBeamTime.push_back(std::make_pair(i, &flash));
156  }
157 
158  for (size_t track_i = 0; track_i < trackVector.size(); track_i++) {
159 
160  recob::Track const& track(trackVector[track_i]);
161 
162  if (track.Length() < fMinTrackLength) continue;
163 
164  //get the begin and end points of this track
165  TVector3 const& pt_begin = track.LocationAtPoint<TVector3>(0);
166  TVector3 const& pt_end = track.LocationAtPoint<TVector3>(track.NumberTrajectoryPoints() - 1);
167  std::vector<float> xyz_begin = {(float)pt_begin.x(), (float)pt_begin.y(), (float)pt_begin.z()};
168  std::vector<float> xyz_end = {(float)pt_end.x(), (float)pt_end.y(), (float)pt_end.z()};
169 
170  //check if this track is outside the drift window, and if it is continue
171  if (!InDriftWindow(pt_begin.x(), pt_end.x(), geom)) continue;
172 
173  cFlashComparison_p.trk_startx = pt_begin.x();
174  cFlashComparison_p.trk_starty = pt_begin.y();
175  cFlashComparison_p.trk_startz = pt_begin.z();
176  cFlashComparison_p.trk_endx = pt_end.x();
177  cFlashComparison_p.trk_endy = pt_end.y();
178  cFlashComparison_p.trk_endz = pt_end.z();
179 
180  //get light hypothesis for track
181  cOpDetVector_hyp = GetMIPHypotheses(track, providers, pvs, opdigip);
182 
183  cFlashComparison_p.hyp_index = track_i;
190  geom);
191 
192  for (auto flash : flashesOnBeamTime) {
193  cOpDetVector_flash = std::vector<float>(geom.NOpDets(), 0);
195  for (size_t c = 0; c <= wireReadoutGeom.MaxOpChannel(); c++) {
196  if (wireReadoutGeom.IsValidOpChannel(c)) {
197  unsigned int OpDet = wireReadoutGeom.OpDetFromOpChannel(c);
198  cOpDetVector_flash[OpDet] += flash.second->PE(c);
199  }
200  }
201  for (size_t o = 0; o < cOpDetVector_flash.size(); o++)
203 
204  cFlashComparison_p.flash_index = flash.first;
205  cFlashComparison_p.flash_totalPE = flash.second->TotalPE();
206  cFlashComparison_p.flash_y = flash.second->YCenter();
207  cFlashComparison_p.flash_sigmay = flash.second->YWidth();
208  cFlashComparison_p.flash_z = flash.second->ZCenter();
209  cFlashComparison_p.flash_sigmaz = flash.second->ZWidth();
210 
212 
213  cTree->Fill();
214  } //end loop over flashes
215 
216  } //end loop over tracks
217 }
218 
219 //this compares the hypothesis to the flash itself.
221  unsigned int const run,
222  unsigned int const event,
223  std::vector<recob::OpFlash> const& flashVector,
224  std::vector<simb::MCParticle> const& mcParticleVector,
225  Providers_t providers,
227  opdet::OpDigiProperties const& opdigip)
228 {
229  auto const& geom = *(providers.get<geo::GeometryCore>());
230  auto const& wireReadoutGeom = art::ServiceHandle<geo::WireReadout const>()->Get();
231 
232  cFlashComparison_p.run = run;
233  cFlashComparison_p.event = event;
234 
235  std::vector<std::pair<unsigned int, const recob::OpFlash*>> flashesOnBeamTime;
236  for (unsigned int i = 0; i < flashVector.size(); i++) {
237  recob::OpFlash const& flash = flashVector[i];
238  if (!flash.OnBeamTime()) continue;
239  flashesOnBeamTime.push_back(std::make_pair(i, &flash));
240  }
241 
242  for (size_t particle_i = 0; particle_i < mcParticleVector.size(); particle_i++) {
243 
244  simb::MCParticle const& particle(mcParticleVector[particle_i]);
245  if (particle.Process().compare("primary") != 0) continue;
246  if (particle.Trajectory().TotalLength() < fMinTrackLength) continue;
247 
248  //get the begin and end points of this track
249  size_t start_i = 0, end_i = particle.NumberTrajectoryPoints() - 1;
250  bool prev_inside = false;
251  for (size_t pt_i = 0; pt_i < particle.NumberTrajectoryPoints(); pt_i++) {
252  bool inside = InDetector(particle.Position(pt_i).Vect(), geom);
253  if (inside && !prev_inside) start_i = pt_i;
254  if (!inside && prev_inside) {
255  end_i = pt_i - 1;
256  break;
257  }
258  prev_inside = inside;
259  }
260  TVector3 const& pt_begin = particle.Position(start_i).Vect();
261  TVector3 const& pt_end = particle.Position(end_i).Vect();
262  std::vector<float> xyz_begin = {(float)pt_begin.x(), (float)pt_begin.y(), (float)pt_begin.z()};
263  std::vector<float> xyz_end = {(float)pt_end.x(), (float)pt_end.y(), (float)pt_end.z()};
264 
265  //check if this track is outside the drift window, and if it is continue
266  if (!InDriftWindow(pt_begin.x(), pt_end.x(), geom)) continue;
267 
268  cFlashComparison_p.trk_startx = pt_begin.x();
269  cFlashComparison_p.trk_starty = pt_begin.y();
270  cFlashComparison_p.trk_startz = pt_begin.z();
271  cFlashComparison_p.trk_endx = pt_end.x();
272  cFlashComparison_p.trk_endy = pt_end.y();
273  cFlashComparison_p.trk_endz = pt_end.z();
274 
275  //get light hypothesis for track
276  cOpDetVector_hyp = GetMIPHypotheses(particle, start_i, end_i, providers, pvs, opdigip);
277 
278  cFlashComparison_p.hyp_index = particle_i;
285  geom);
286 
287  for (auto flash : flashesOnBeamTime) {
288  cOpDetVector_flash = std::vector<float>(geom.NOpDets(), 0);
290  for (size_t c = 0; c <= wireReadoutGeom.MaxOpChannel(); c++) {
291  if (wireReadoutGeom.IsValidOpChannel(c)) {
292  unsigned int OpDet = wireReadoutGeom.OpDetFromOpChannel(c);
293  cOpDetVector_flash[OpDet] += flash.second->PE(c);
294  }
295  }
296  for (size_t o = 0; o < cOpDetVector_flash.size(); o++)
298  cFlashComparison_p.flash_index = flash.first;
299  cFlashComparison_p.flash_totalPE = flash.second->TotalPE();
300  cFlashComparison_p.flash_y = flash.second->YCenter();
301  cFlashComparison_p.flash_sigmay = flash.second->YWidth();
302  cFlashComparison_p.flash_z = flash.second->ZCenter();
303  cFlashComparison_p.flash_sigmaz = flash.second->ZWidth();
304 
306 
307  for (size_t i = 0; i < cOpDetVector_flash.size(); i++)
308  cOpDetHist_flash->SetBinContent(i + 1, cOpDetVector_flash[i]);
309 
310  for (size_t i = 0; i < cOpDetVector_hyp.size(); i++)
311  cOpDetHist_hyp->SetBinContent(i + 1, cOpDetVector_hyp[i]);
312 
313  for (size_t i = 0; i < cOpDetVector_flash.size(); i++) {
314  std::cout << "Flash/Hyp " << i << " : " << cOpDetHist_flash->GetBinContent(i + 1) << " "
315  << cOpDetHist_hyp->GetBinContent(i + 1) << std::endl;
316  }
317 
318  cTree->Fill();
319  } //end loop over flashes
320 
321  } //end loop over tracks
322 }
323 
325  std::vector<float> const& opdetVector,
326  float& sum,
327  float& y,
328  float& sigmay,
329  float& z,
330  float& sigmaz,
331  geo::GeometryCore const& geom)
332 {
333  y = 0;
334  sigmay = 0;
335  z = 0;
336  sigmaz = 0;
337  sum = 0;
338  for (unsigned int opdet = 0; opdet < opdetVector.size(); opdet++) {
339  sum += opdetVector[opdet];
340  auto const xyz = geom.Cryostat().OpDet(opdet).GetCenter();
341  y += opdetVector[opdet] * xyz.Y();
342  z += opdetVector[opdet] * xyz.Z();
343  }
344 
345  y /= sum;
346  z /= sum;
347 
348  for (unsigned int opdet = 0; opdet < opdetVector.size(); opdet++) {
349  auto const xyz = geom.Cryostat().OpDet(opdet).GetCenter();
350  sigmay += (opdetVector[opdet] * xyz.Y() - y) * (opdetVector[opdet] * xyz.Y() - y);
351  sigmaz += (opdetVector[opdet] * xyz.Z() - y) * (opdetVector[opdet] * xyz.Z() - y);
352  }
353 
354  sigmay = std::sqrt(sigmay) / sum;
355  sigmaz = std::sqrt(sigmaz) / sum;
356 }
357 
359  geo::GeometryCore const& geom)
360 {
361  auto const& tpc = geom.TPC({0, 0});
362  if (pt.x() < 0 || pt.x() > 2 * tpc.HalfWidth()) return false;
363  if (std::abs(pt.y()) > tpc.HalfHeight()) return false;
364  if (pt.z() < 0 || pt.z() > tpc.Length()) return false;
365  return true;
366 }
367 
369  double end_x,
370  geo::GeometryCore const& geom)
371 {
372  if (start_x < 0. || end_x < 0.) return false;
373  auto const& tpc = geom.TPC({0, 0});
374  if (start_x > 2 * tpc.HalfWidth() || end_x > 2 * tpc.HalfWidth()) return false;
375  return true;
376 }
377 
379  TVector3 const& pt1,
380  TVector3 const& pt2,
381  std::vector<float>& lightHypothesis,
382  float& totalHypothesisPE,
383  geo::GeometryCore const& geom,
385  float const& PromptMIPScintYield,
386  float XOffset)
387 {
388  double xyz_segment[3];
389  xyz_segment[0] = 0.5 * (pt2.x() + pt1.x()) + XOffset;
390  xyz_segment[1] = 0.5 * (pt2.y() + pt1.y());
391  xyz_segment[2] = 0.5 * (pt2.z() + pt1.z());
392 
393  //get the visibility vector
394  auto const& PointVisibility = pvs.GetAllVisibilities(xyz_segment);
395 
396  //check pointer, as it may be null if given a y/z outside some range
397  if (!PointVisibility) return;
398 
399  //get the amount of light
400  float LightAmount = PromptMIPScintYield * (pt2 - pt1).Mag();
401 
402  for (size_t opdet_i = 0; opdet_i < pvs.NOpChannels(); opdet_i++) {
403  lightHypothesis[opdet_i] += PointVisibility[opdet_i] * LightAmount;
404  totalHypothesisPE += PointVisibility[opdet_i] * LightAmount;
405 
406  //apply saturation limit
407  if (lightHypothesis[opdet_i] > fOpDetSaturation) {
408  totalHypothesisPE -= (lightHypothesis[opdet_i] - fOpDetSaturation);
409  lightHypothesis[opdet_i] = fOpDetSaturation;
410  }
411  }
412 
413 } //end AddLightFromSegment
414 
416  std::vector<float>& lightHypothesis,
417  float const& totalHypothesisPE,
418  geo::GeometryCore const& geom)
419 {
420  for (size_t opdet_i = 0; opdet_i < geom.NOpDets(); opdet_i++)
421  lightHypothesis[opdet_i] /= totalHypothesisPE;
422 }
423 
424 // Get a hypothesis for the light collected for a track
426  recob::Track const& track,
427  Providers_t providers,
429  opdet::OpDigiProperties const& opdigip,
430  float XOffset)
431 {
432  auto const& geom = *(providers.get<geo::GeometryCore>());
433  auto const& larp = *(providers.get<detinfo::LArProperties>());
434  std::vector<float> lightHypothesis(geom.NOpDets(), 0);
435  float totalHypothesisPE = 0;
436  const float PromptMIPScintYield =
437  larp.ScintYield() * larp.ScintYieldRatio() * opdigip.QE() * fMIPdQdx;
438 
439  //get QE from ubChannelConfig, which gives per tube, so goes in AddLightFromSegment
440  //VisibleEnergySeparation(step);
441 
442  for (size_t pt = 1; pt < track.NumberTrajectoryPoints(); pt++)
443  AddLightFromSegment(track.LocationAtPoint<TVector3>(pt - 1),
444  track.LocationAtPoint<TVector3>(pt),
445  lightHypothesis,
446  totalHypothesisPE,
447  geom,
448  pvs,
449  PromptMIPScintYield,
450  XOffset);
451 
452  if (fNormalizeHypothesisToFlash && totalHypothesisPE > std::numeric_limits<float>::epsilon())
453  NormalizeLightHypothesis(lightHypothesis, totalHypothesisPE, geom);
454 
455  return lightHypothesis;
456 
457 } //end GetMIPHypotheses
458 
459 // Get a hypothesis for the light collected for a particle trajectory
461  simb::MCParticle const& particle,
462  size_t start_i,
463  size_t end_i,
464  Providers_t providers,
466  opdet::OpDigiProperties const& opdigip,
467  float XOffset)
468 {
469  auto const& geom = *(providers.get<geo::GeometryCore>());
470  auto const& larp = *(providers.get<detinfo::LArProperties>());
471  std::vector<float> lightHypothesis(geom.NOpDets(), 0);
472  float totalHypothesisPE = 0;
473  const float PromptMIPScintYield =
474  larp.ScintYield() * larp.ScintYieldRatio() * opdigip.QE() * fMIPdQdx;
475 
476  for (size_t pt = start_i + 1; pt <= end_i; pt++)
477  AddLightFromSegment(particle.Position(pt - 1).Vect(),
478  particle.Position(pt).Vect(),
479  lightHypothesis,
480  totalHypothesisPE,
481  geom,
482  pvs,
483  PromptMIPScintYield,
484  XOffset);
485 
486  if (fNormalizeHypothesisToFlash && totalHypothesisPE > std::numeric_limits<float>::epsilon())
487  NormalizeLightHypothesis(lightHypothesis, totalHypothesisPE, geom);
488 
489  return lightHypothesis;
490 
491 } //end GetMIPHypotheses
492 
493 //---------------------------------------
494 // Check whether a hypothesis can be accomodated in a flash
495 // Flashes fail if 1 bin is far in excess of the observed signal
496 // or if the whole flash intensity is much too large for the hypothesis.
497 // MIP dEdx is assumed for now. Accounting for real dQdx will
498 // improve performance of this algorithm.
499 //---------------------------------------
502  std::vector<float> const& lightHypothesis,
503  const recob::OpFlash* flashPointer,
504  geo::GeometryCore const& geom,
505  geo::WireReadoutGeom const& wireReadoutGeom)
506 {
507  float hypothesis_integral = 0;
508  float flash_integral = 0;
509  unsigned int cumulativeChannels = 0;
510 
511  std::vector<double> PEbyOpDet(geom.NOpDets(), 0);
512  for (unsigned int c = 0; c <= wireReadoutGeom.MaxOpChannel(); c++) {
513  if (wireReadoutGeom.IsValidOpChannel(c)) {
514  unsigned int o = wireReadoutGeom.OpDetFromOpChannel(c);
515  PEbyOpDet[o] += flashPointer->PE(c);
516  }
517  }
518 
519  float hypothesis_scale = 1.;
520  if (fNormalizeHypothesisToFlash) hypothesis_scale = flashPointer->TotalPE();
521 
522  for (size_t pmt_i = 0; pmt_i < lightHypothesis.size(); pmt_i++) {
523 
524  flash_integral += PEbyOpDet[pmt_i];
525 
526  if (lightHypothesis[pmt_i] < std::numeric_limits<float>::epsilon()) continue;
527  hypothesis_integral += lightHypothesis[pmt_i] * hypothesis_scale;
528 
529  if (PEbyOpDet[pmt_i] < fMinOpHitPE) continue;
530 
531  float diff_scaled = (lightHypothesis[pmt_i] * hypothesis_scale - PEbyOpDet[pmt_i]) /
532  std::sqrt(lightHypothesis[pmt_i] * hypothesis_scale);
533 
534  if (diff_scaled > fSingleChannelCut) return CompatibilityResultType::kSingleChannelCut;
535 
536  if (diff_scaled > fCumulativeChannelThreshold) cumulativeChannels++;
537  if (cumulativeChannels >= fCumulativeChannelCut)
538  return CompatibilityResultType::kCumulativeChannelCut;
539  }
540 
541  if ((hypothesis_integral - flash_integral) / std::sqrt(hypothesis_integral) > fIntegralCut)
542  return CompatibilityResultType::kIntegralCut;
543 
544  return CompatibilityResultType::kCompatible;
545 }
546 
547 float cosmic::BeamFlashTrackMatchTaggerAlg::CalculateChi2(std::vector<float> const& light_flash,
548  std::vector<float> const& light_track)
549 {
550  float chi2 = 0;
551  for (size_t pmt_i = 0; pmt_i < light_flash.size(); pmt_i++) {
552 
553  if (light_flash[pmt_i] < fMinOpHitPE) continue;
554 
555  float err2 = 1;
556  if (light_track[pmt_i] > 1) err2 = light_track[pmt_i];
557 
558  chi2 +=
559  (light_flash[pmt_i] - light_track[pmt_i]) * (light_flash[pmt_i] - light_track[pmt_i]) / err2;
560  }
561 
562  return chi2;
563 }
564 
566  std::ostream* output)
567 {
568  *output << "----------------------------------------------" << std::endl;
569  *output << "Track properties: ";
570  *output << "\n\tLength=" << track.Length();
571 
572  auto const& pt_begin = track.LocationAtPoint(0);
573  *output << "\n\tBegin Location (x,y,z)=(" << pt_begin.x() << "," << pt_begin.y() << ","
574  << pt_begin.z() << ")";
575 
576  auto const& pt_end = track.LocationAtPoint(track.NumberTrajectoryPoints() - 1);
577  *output << "\n\tEnd Location (x,y,z)=(" << pt_end.x() << "," << pt_end.y() << "," << pt_end.z()
578  << ")";
579 
580  *output << "\n\tTrajectoryPoints=" << track.NumberTrajectoryPoints();
581  *output << std::endl;
582  *output << "----------------------------------------------" << std::endl;
583 }
584 
586  std::ostream* output)
587 {
588  *output << "----------------------------------------------" << std::endl;
589  *output << "Flash properties: ";
590 
591  *output << "\n\tTime=" << flash.Time();
592  *output << "\n\tOnBeamTime=" << flash.OnBeamTime();
593  *output << "\n\ty position (center,width)=(" << flash.YCenter() << "," << flash.YWidth() << ")";
594  *output << "\n\tz position (center,width)=(" << flash.ZCenter() << "," << flash.ZWidth() << ")";
595  *output << "\n\tTotal PE=" << flash.TotalPE();
596 
597  *output << std::endl;
598  *output << "----------------------------------------------" << std::endl;
599 }
600 
602  std::vector<float> const& lightHypothesis,
603  const recob::OpFlash* flashPointer,
604  geo::GeometryCore const& geom,
605  geo::WireReadoutGeom const& wireReadoutGeom,
607  std::ostream* output)
608 {
609  *output << "----------------------------------------------" << std::endl;
610  *output << "Hypothesis-flash comparison: ";
611 
612  float hypothesis_integral = 0;
613  float flash_integral = 0;
614 
615  float hypothesis_scale = 1.;
616  if (fNormalizeHypothesisToFlash) hypothesis_scale = flashPointer->TotalPE();
617 
618  std::vector<double> PEbyOpDet(geom.NOpDets(), 0);
619  for (unsigned int c = 0; c <= wireReadoutGeom.MaxOpChannel(); c++) {
620  if (wireReadoutGeom.IsValidOpChannel(c)) {
621  unsigned int o = wireReadoutGeom.OpDetFromOpChannel(c);
622  PEbyOpDet[o] += flashPointer->PE(c);
623  }
624  }
625 
626  for (size_t pmt_i = 0; pmt_i < lightHypothesis.size(); pmt_i++) {
627 
628  flash_integral += PEbyOpDet[pmt_i];
629 
630  *output << "\n\t pmt_i=" << pmt_i << ", (hypothesis,flash)=("
631  << lightHypothesis[pmt_i] * hypothesis_scale << "," << PEbyOpDet[pmt_i] << ")";
632 
633  if (lightHypothesis[pmt_i] < std::numeric_limits<float>::epsilon()) continue;
634 
635  *output << " difference="
636  << (lightHypothesis[pmt_i] * hypothesis_scale - PEbyOpDet[pmt_i]) /
637  std::sqrt(lightHypothesis[pmt_i] * hypothesis_scale);
638 
639  hypothesis_integral += lightHypothesis[pmt_i] * hypothesis_scale;
640  }
641 
642  *output << "\n\t TOTAL (hypothesis,flash)=(" << hypothesis_integral << "," << flash_integral
643  << ")"
644  << " difference="
645  << (hypothesis_integral - flash_integral) / std::sqrt(hypothesis_integral);
646 
647  *output << std::endl;
648  *output << "End result=" << result << std::endl;
649  *output << "----------------------------------------------" << std::endl;
650 }
unsigned int NumberTrajectoryPoints() const
Definition: MCParticle.h:219
const TLorentzVector & Position(const int i=0) const
Definition: MCParticle.h:220
CompatibilityResultType CheckCompatibility(std::vector< float > const &lightHypothesis, const recob::OpFlash *flashPointer, geo::GeometryCore const &geom, geo::WireReadoutGeom const &wireReadoutGeom)
BeamFlashTrackMatchTaggerAlg(fhicl::ParameterSet const &p)
void PrintHypothesisFlashComparison(std::vector< float > const &, const recob::OpFlash *, geo::GeometryCore const &geom, geo::WireReadoutGeom const &wireReadoutGeom, CompatibilityResultType, std::ostream *output=&std::cout)
OpDetGeo const & OpDet(unsigned int iopdet) const
Return the iopdet&#39;th optical detector in the cryostat.
Definition: CryostatGeo.cxx:83
const simb::MCTrajectory & Trajectory() const
Definition: MCParticle.h:254
bool IsValidOpChannel(int opChannel) const
Is this a valid OpChannel number?
enum anab::cosmic_tag_id CosmicTagID_t
Float_t y
Definition: compare.C:6
Point_t const & GetCenter() const
Definition: OpDetGeo.h:71
Point_t const & LocationAtPoint(size_t i) const
Access to track position at different points.
Definition: Track.h:160
Double_t z
Definition: plot.C:276
void RunHypothesisComparison(unsigned int const, unsigned int const, std::vector< recob::OpFlash > const &, std::vector< recob::Track > const &, Providers_t, phot::PhotonVisibilityService const &, opdet::OpDigiProperties const &)
float CalculateChi2(std::vector< float > const &, std::vector< float > const &)
size_t NumberTrajectoryPoints() const
Various functions related to the presence and the number of (valid) points.
Definition: Track.h:136
constexpr auto abs(T v)
Returns the absolute value of the argument.
virtual unsigned int OpDetFromOpChannel(unsigned int opChannel) const
Returns the optical detector the specified optical channel belongs.
Provider const * get() const
Returns the provider with the specified type.
Definition: ProviderPack.h:182
cout<< "Opened file "<< fin<< " ixs= "<< ixs<< endl;if(ixs==0) hhh=(TH1F *) fff-> Get("h1")
Definition: AddMC.C:8
double PE(unsigned int i) const
Definition: OpFlash.h:134
std::string Process() const
Definition: MCParticle.h:216
void RunCompatibilityCheck(std::vector< recob::OpFlash > const &, std::vector< recob::Track > const &, std::vector< anab::CosmicTag > &, std::vector< size_t > &, Providers_t, phot::PhotonVisibilityService const &, opdet::OpDigiProperties const &)
double ZCenter() const
Definition: OpFlash.h:162
Access the description of the physical detector geometry.
std::vector< float > GetMIPHypotheses(recob::Track const &track, Providers_t providers, phot::PhotonVisibilityService const &pvs, opdet::OpDigiProperties const &, float XOffset=0)
double Time() const
Definition: OpFlash.h:118
double QE() const noexcept
Returns quantum efficiency.
double Length(size_t p=0) const
Access to various track properties.
Definition: Track.h:207
TText * pt2
Definition: plot.C:64
void NormalizeLightHypothesis(std::vector< float > &lightHypothesis, float const &totalHypothesisPE, geo::GeometryCore const &geom)
int OnBeamTime() const
Definition: OpFlash.h:186
void AddLightFromSegment(TVector3 const &pt1, TVector3 const &pt2, std::vector< float > &lightHypothesis, float &totalHypothesisPE, geo::GeometryCore const &geom, phot::PhotonVisibilityService const &pvs, float const &PromptMIPScintYield, float XOffset)
Interface for a class providing readout channel mapping to geometry.
void PrintTrackProperties(recob::Track const &, std::ostream *output=&std::cout)
TMarker * pt
Definition: egs.C:25
Description of the physical geometry of one entire detector.
Definition: GeometryCore.h:91
MappedCounts_t GetAllVisibilities(Point const &p, bool wantReflected=false) const
void SetHypothesisComparisonTree(TTree *, TH1F *, TH1F *)
CryostatGeo const & Cryostat(CryostatID const &cryoid=details::cryostat_zero) const
Returns the specified cryostat.
unsigned int NOpDets() const
Number of OpDets in the whole detector.
Encapsulate the geometry of an optical detector.
double YWidth() const
Definition: OpFlash.h:158
decltype(auto) get(T &&obj)
ADL-aware version of std::to_string.
Definition: StdUtils.h:120
double TotalLength() const
unsigned int MaxOpChannel() const
Largest optical channel number.
Container for a list of pointers to providers.
Definition: ProviderPack.h:111
bool InDriftWindow(double, double, geo::GeometryCore const &)
void FillFlashProperties(std::vector< float > const &opdetVector, float &, float &, float &, float &, float &, geo::GeometryCore const &geom)
double TotalPE() const
Definition: OpFlash.cxx:94
TText * pt1
Definition: plot.C:61
bool InDetector(TVector3 const &, geo::GeometryCore const &)
TPCGeo const & TPC(TPCID const &tpcid=details::tpc_zero) const
Returns the specified TPC.
Definition: GeometryCore.h:448
Double_t sum
Definition: plot.C:31
double YCenter() const
Definition: OpFlash.h:154
Float_t track
Definition: plot.C:35
art framework interface to geometry description
Track from a non-cascading particle.A recob::Track consists of a recob::TrackTrajectory, plus additional members relevant for a "fitted" track:
Definition: Track.h:49
double ZWidth() const
Definition: OpFlash.h:166
Event finding and building.
void PrintFlashProperties(recob::OpFlash const &, std::ostream *output=&std::cout)