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"))
55 cOpDetHist_flash->SetNameTitle(
"opdet_hist_flash",
"Optical Detector Occupancy, Flash");
58 cOpDetHist_hyp->SetNameTitle(
"opdet_hist_hyp",
"Optical Detector Occupancy, Hypothesis");
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,
79 std::vector<const recob::OpFlash*> flashesOnBeamTime;
80 for (
auto const& flash : flashVector) {
81 if (!flash.OnBeamTime())
continue;
82 flashesOnBeamTime.push_back(&flash);
86 assnTrackTagVector.resize(trackVector.size(), std::numeric_limits<size_t>::max());
87 cosmicTagVector.reserve(trackVector.size());
89 for (
size_t track_i = 0; track_i < trackVector.size(); track_i++) {
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()};
105 assnTrackTagVector[track_i] = cosmicTagVector.size() - 1;
111 std::vector<float> lightHypothesis =
GetMIPHypotheses(track, providers, pvs, opdigip);
114 bool compatible =
false;
118 if (result == CompatibilityResultType::kCompatible) compatible =
true;
123 lightHypothesis, flashPointer, geom, wireReadoutGeom, result);
128 float cosmicScore = 1.;
129 if (compatible) cosmicScore = 0.;
131 assnTrackTagVector[track_i] = cosmicTagVector.size() - 1;
137 unsigned int const run,
138 unsigned int const event,
139 std::vector<recob::OpFlash>
const& flashVector,
140 std::vector<recob::Track>
const& trackVector,
151 std::vector<std::pair<unsigned int, const recob::OpFlash*>> flashesOnBeamTime;
152 for (
unsigned int i = 0; i < flashVector.size(); i++) {
155 flashesOnBeamTime.push_back(std::make_pair(i, &flash));
158 for (
size_t track_i = 0; track_i < trackVector.size(); track_i++) {
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()};
171 if (!
InDriftWindow(pt_begin.x(), pt_end.x(), geom))
continue;
192 for (
auto flash : flashesOnBeamTime) {
195 for (
size_t c = 0; c <= wireReadoutGeom.MaxOpChannel(); c++) {
196 if (wireReadoutGeom.IsValidOpChannel(c)) {
197 unsigned int OpDet = wireReadoutGeom.OpDetFromOpChannel(c);
221 unsigned int const run,
222 unsigned int const event,
223 std::vector<recob::OpFlash>
const& flashVector,
224 std::vector<simb::MCParticle>
const& mcParticleVector,
235 std::vector<std::pair<unsigned int, const recob::OpFlash*>> flashesOnBeamTime;
236 for (
unsigned int i = 0; i < flashVector.size(); i++) {
239 flashesOnBeamTime.push_back(std::make_pair(i, &flash));
242 for (
size_t particle_i = 0; particle_i < mcParticleVector.size(); particle_i++) {
245 if (particle.
Process().compare(
"primary") != 0)
continue;
250 bool prev_inside =
false;
253 if (inside && !prev_inside) start_i = pt_i;
254 if (!inside && prev_inside) {
258 prev_inside = inside;
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()};
266 if (!
InDriftWindow(pt_begin.x(), pt_end.x(), geom))
continue;
287 for (
auto flash : flashesOnBeamTime) {
290 for (
size_t c = 0; c <= wireReadoutGeom.MaxOpChannel(); c++) {
291 if (wireReadoutGeom.IsValidOpChannel(c)) {
292 unsigned int OpDet = wireReadoutGeom.OpDetFromOpChannel(c);
314 std::cout <<
"Flash/Hyp " << i <<
" : " <<
cOpDetHist_flash->GetBinContent(i + 1) <<
" " 325 std::vector<float>
const& opdetVector,
339 sum += opdetVector[
opdet];
341 y += opdetVector[
opdet] * xyz.Y();
342 z += opdetVector[
opdet] * xyz.Z();
350 sigmay += (opdetVector[
opdet] * xyz.Y() -
y) * (opdetVector[
opdet] * xyz.Y() -
y);
351 sigmaz += (opdetVector[
opdet] * xyz.Z() -
y) * (opdetVector[
opdet] * xyz.Z() -
y);
354 sigmay = std::sqrt(sigmay) /
sum;
355 sigmaz = std::sqrt(sigmaz) /
sum;
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;
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;
381 std::vector<float>& lightHypothesis,
382 float& totalHypothesisPE,
385 float const& PromptMIPScintYield,
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());
397 if (!PointVisibility)
return;
400 float LightAmount = PromptMIPScintYield * (pt2 -
pt1).Mag();
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;
416 std::vector<float>& lightHypothesis,
417 float const& totalHypothesisPE,
420 for (
size_t opdet_i = 0; opdet_i < geom.
NOpDets(); opdet_i++)
421 lightHypothesis[opdet_i] /= totalHypothesisPE;
434 std::vector<float> lightHypothesis(geom.NOpDets(), 0);
435 float totalHypothesisPE = 0;
436 const float PromptMIPScintYield =
437 larp.ScintYield() * larp.ScintYieldRatio() * opdigip.
QE() *
fMIPdQdx;
455 return lightHypothesis;
471 std::vector<float> lightHypothesis(geom.NOpDets(), 0);
472 float totalHypothesisPE = 0;
473 const float PromptMIPScintYield =
474 larp.ScintYield() * larp.ScintYieldRatio() * opdigip.
QE() *
fMIPdQdx;
476 for (
size_t pt = start_i + 1;
pt <= end_i;
pt++)
489 return lightHypothesis;
502 std::vector<float>
const& lightHypothesis,
507 float hypothesis_integral = 0;
508 float flash_integral = 0;
509 unsigned int cumulativeChannels = 0;
511 std::vector<double> PEbyOpDet(geom.
NOpDets(), 0);
512 for (
unsigned int c = 0; c <= wireReadoutGeom.
MaxOpChannel(); c++) {
515 PEbyOpDet[o] += flashPointer->
PE(c);
519 float hypothesis_scale = 1.;
522 for (
size_t pmt_i = 0; pmt_i < lightHypothesis.size(); pmt_i++) {
524 flash_integral += PEbyOpDet[pmt_i];
526 if (lightHypothesis[pmt_i] < std::numeric_limits<float>::epsilon())
continue;
527 hypothesis_integral += lightHypothesis[pmt_i] * hypothesis_scale;
531 float diff_scaled = (lightHypothesis[pmt_i] * hypothesis_scale - PEbyOpDet[pmt_i]) /
532 std::sqrt(lightHypothesis[pmt_i] * hypothesis_scale);
534 if (diff_scaled >
fSingleChannelCut)
return CompatibilityResultType::kSingleChannelCut;
538 return CompatibilityResultType::kCumulativeChannelCut;
541 if ((hypothesis_integral - flash_integral) / std::sqrt(hypothesis_integral) >
fIntegralCut)
542 return CompatibilityResultType::kIntegralCut;
544 return CompatibilityResultType::kCompatible;
548 std::vector<float>
const& light_track)
551 for (
size_t pmt_i = 0; pmt_i < light_flash.size(); pmt_i++) {
556 if (light_track[pmt_i] > 1) err2 = light_track[pmt_i];
559 (light_flash[pmt_i] - light_track[pmt_i]) * (light_flash[pmt_i] - light_track[pmt_i]) / err2;
566 std::ostream* output)
568 *output <<
"----------------------------------------------" << std::endl;
569 *output <<
"Track properties: ";
570 *output <<
"\n\tLength=" << track.
Length();
573 *output <<
"\n\tBegin Location (x,y,z)=(" << pt_begin.x() <<
"," << pt_begin.y() <<
"," 574 << pt_begin.z() <<
")";
577 *output <<
"\n\tEnd Location (x,y,z)=(" << pt_end.x() <<
"," << pt_end.y() <<
"," << pt_end.z()
581 *output << std::endl;
582 *output <<
"----------------------------------------------" << std::endl;
586 std::ostream* output)
588 *output <<
"----------------------------------------------" << std::endl;
589 *output <<
"Flash properties: ";
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();
597 *output << std::endl;
598 *output <<
"----------------------------------------------" << std::endl;
602 std::vector<float>
const& lightHypothesis,
607 std::ostream* output)
609 *output <<
"----------------------------------------------" << std::endl;
610 *output <<
"Hypothesis-flash comparison: ";
612 float hypothesis_integral = 0;
613 float flash_integral = 0;
615 float hypothesis_scale = 1.;
618 std::vector<double> PEbyOpDet(geom.
NOpDets(), 0);
619 for (
unsigned int c = 0; c <= wireReadoutGeom.
MaxOpChannel(); c++) {
622 PEbyOpDet[o] += flashPointer->
PE(c);
626 for (
size_t pmt_i = 0; pmt_i < lightHypothesis.size(); pmt_i++) {
628 flash_integral += PEbyOpDet[pmt_i];
630 *output <<
"\n\t pmt_i=" << pmt_i <<
", (hypothesis,flash)=(" 631 << lightHypothesis[pmt_i] * hypothesis_scale <<
"," << PEbyOpDet[pmt_i] <<
")";
633 if (lightHypothesis[pmt_i] < std::numeric_limits<float>::epsilon())
continue;
635 *output <<
" difference=" 636 << (lightHypothesis[pmt_i] * hypothesis_scale - PEbyOpDet[pmt_i]) /
637 std::sqrt(lightHypothesis[pmt_i] * hypothesis_scale);
639 hypothesis_integral += lightHypothesis[pmt_i] * hypothesis_scale;
642 *output <<
"\n\t TOTAL (hypothesis,flash)=(" << hypothesis_integral <<
"," << flash_integral
645 << (hypothesis_integral - flash_integral) / std::sqrt(hypothesis_integral);
647 *output << std::endl;
648 *output <<
"End result=" << result << std::endl;
649 *output <<
"----------------------------------------------" << std::endl;
unsigned int NumberTrajectoryPoints() const
const TLorentzVector & Position(const int i=0) const
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'th optical detector in the cryostat.
float fCumulativeChannelThreshold
const simb::MCTrajectory & Trajectory() const
bool IsValidOpChannel(int opChannel) const
Is this a valid OpChannel number?
enum anab::cosmic_tag_id CosmicTagID_t
Point_t const & GetCenter() const
Point_t const & LocationAtPoint(size_t i) const
Access to track position at different points.
bool fNormalizeHypothesisToFlash
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.
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.
cout<< "Opened file "<< fin<< " ixs= "<< ixs<< endl;if(ixs==0) hhh=(TH1F *) fff-> Get("h1")
double PE(unsigned int i) const
std::string Process() const
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 &)
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)
const anab::CosmicTagID_t COSMIC_TYPE_OUTSIDEDRIFT
double QE() const noexcept
Returns quantum efficiency.
std::vector< float > cOpDetVector_flash
double Length(size_t p=0) const
Access to various track properties.
void NormalizeLightHypothesis(std::vector< float > &lightHypothesis, float const &totalHypothesisPE, geo::GeometryCore const &geom)
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.
FlashComparisonProperties_t cFlashComparison_p
bool fMakeOutsideDriftTags
void PrintTrackProperties(recob::Track const &, std::ostream *output=&std::cout)
unsigned int flash_nOpDet
Description of the physical geometry of one entire detector.
std::string leaf_structure
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.
decltype(auto) get(T &&obj)
ADL-aware version of std::to_string.
unsigned int fCumulativeChannelCut
double TotalLength() const
std::vector< float > cOpDetVector_hyp
unsigned int MaxOpChannel() const
Largest optical channel number.
Container for a list of pointers to providers.
bool InDriftWindow(double, double, geo::GeometryCore const &)
void FillFlashProperties(std::vector< float > const &opdetVector, float &, float &, float &, float &, float &, geo::GeometryCore const &geom)
const anab::CosmicTagID_t COSMIC_TYPE_FLASHMATCH
bool InDetector(TVector3 const &, geo::GeometryCore const &)
size_t NOpChannels() const
TPCGeo const & TPC(TPCID const &tpcid=details::tpc_zero) const
Returns the specified TPC.
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:
Event finding and building.
void PrintFlashProperties(recob::OpFlash const &, std::ostream *output=&std::cout)