File indexing completed on 2025-08-06 08:18:10
0001 #ifndef TRACKBASE_RESIDUALOUTLIERFINDER_H
0002 #define TRACKBASE_RESIDUALOUTLIERFINDER_H
0003
0004 #include <TFile.h>
0005 #include <TH2.h>
0006 #include <TNtuple.h>
0007 #include <phool/phool.h>
0008 #include <Acts/Definitions/Units.hpp>
0009 #include <Acts/EventData/Measurement.hpp>
0010 #include <Acts/EventData/MeasurementHelpers.hpp>
0011 #include <Acts/EventData/MultiTrajectory.hpp>
0012 #include <Acts/EventData/VectorMultiTrajectory.hpp>
0013
0014 #include <Acts/EventData/MultiTrajectoryHelpers.hpp>
0015 struct ResidualOutlierFinder
0016 {
0017 ActsGeometry* m_tGeometry = nullptr;
0018 int verbosity = 0;
0019 std::map<long unsigned int, float> chi2Cuts;
0020 std::string outfilename = "OutlierFinder.root";
0021 TH2 *hChi2 = new TH2F("h_layerChi2", ";layer;#chi^{2}", 60, 0, 60, 1000, 0, 500);
0022
0023 TH2 *hDistance = new TH2F("h_layerDistance", ";layer;distance", 60, 0, 60, 1000, 0, 500);
0024
0025 TNtuple *tree = new TNtuple("ntp_outlierfinder", "tree with predicted states and cluster info",
0026 "sphenixlayer:layer:volume:distance:chi2:predgx:predgy:predgz:predlx:predlz:clusgx:clusgy:clusgz:cluslx:cluslz");
0027 void outfileName(const std::string& name)
0028 {
0029 outfilename = name;
0030 }
0031 void Write()
0032 {
0033 TFile *file = new TFile(outfilename.c_str(), "recreate");
0034
0035 file->cd();
0036 hChi2->Write();
0037 hDistance->Write();
0038 tree->Write();
0039 file->ls();
0040 }
0041 bool operator()(Acts::MultiTrajectory<Acts::VectorMultiTrajectory>::ConstTrackStateProxy state) const
0042 {
0043
0044 if (!state.hasCalibrated() || !state.hasPredicted())
0045 {
0046 return false;
0047 }
0048
0049 auto sourceLink = state.getUncalibratedSourceLink().template get<ActsSourceLink>();
0050 const auto& cluskey = sourceLink.cluskey();
0051
0052 const auto predicted = state.predicted();
0053 const auto predictedCovariance = state.predictedCovariance();
0054 float chi2 = std::numeric_limits<float>::max();
0055
0056 auto fullCalibrated = state
0057 .template calibrated<Acts::MultiTrajectoryTraits::MeasurementSizeMax>()
0058 .data();
0059 auto fullCalibratedCovariance = state
0060 .template calibratedCovariance<Acts::MultiTrajectoryTraits::MeasurementSizeMax>()
0061 .data();
0062
0063 chi2 = Acts::visit_measurement(state.calibratedSize(), [&](auto N) -> double
0064 {
0065 constexpr size_t kMeasurementSize = decltype(N)::value;
0066 typename Acts::TrackStateTraits<kMeasurementSize, true>::Measurement calibrated{
0067 fullCalibrated};
0068
0069 typename Acts::TrackStateTraits<kMeasurementSize, true>::MeasurementCovariance
0070 calibratedCovariance{fullCalibratedCovariance};
0071
0072 using ParametersVector = Acts::ActsVector<kMeasurementSize>;
0073 const auto H = state.projector().template topLeftCorner<kMeasurementSize, Acts::eBoundSize>().eval();
0074 ParametersVector res;
0075 res = calibrated - H * predicted;
0076 chi2 = (res.transpose() * ((calibratedCovariance + H * predictedCovariance * H.transpose())).inverse() * res).eval()(0, 0);
0077
0078 return chi2; });
0079
0080 float distance = Acts::visit_measurement(state.calibratedSize(), [&](auto N)
0081 {
0082 constexpr size_t kMeasurementSize = decltype(N)::value;
0083 auto residuals =
0084 state.template calibrated<kMeasurementSize>() -
0085 state.projector()
0086 .template topLeftCorner<kMeasurementSize, Acts::eBoundSize>() *
0087 state.predicted();
0088 auto cdistance = residuals.norm();
0089 return cdistance; });
0090
0091 if (verbosity > 2)
0092 {
0093 std::cout << "Measurement has distance, chi2 "
0094 << distance << ", " << chi2
0095 << std::endl;
0096 }
0097
0098 auto volume = state.referenceSurface().geometryId().volume();
0099 auto layer = state.referenceSurface().geometryId().layer();
0100
0101 bool outlier = false;
0102 float chi2cut = chi2Cuts.find(volume)->second;
0103 int sphenixlayer = TrkrDefs::getLayer(cluskey);
0104 hChi2->Fill(sphenixlayer, chi2);
0105 hDistance->Fill(sphenixlayer, distance);
0106 if(!m_tGeometry)
0107 {
0108 std::cout << PHWHERE << "no geometry set in residual outlier finder" << std::endl;
0109 exit(1);
0110 }
0111 Acts::FreeVector freeParams =
0112 Acts::transformBoundToFreeParameters(state.referenceSurface(),
0113 m_tGeometry->geometry().getGeoContext(),
0114 predicted);
0115 Acts::Vector2 local(fullCalibrated[Acts::eBoundLoc0], fullCalibrated[Acts::eBoundLoc1]);
0116 Acts::Vector3 global = state.referenceSurface().localToGlobal(
0117 m_tGeometry->geometry().getGeoContext(), local,
0118 Acts::Vector3(1, 1, 1));
0119 float data[] = {
0120 (float) sphenixlayer, (float) layer, (float) volume, distance, chi2,
0121 (float) freeParams[Acts::eFreePos0], (float) freeParams[Acts::eFreePos1], (float) freeParams[Acts::eFreePos2],
0122 (float) predicted[Acts::eBoundLoc0], (float) predicted[Acts::eBoundLoc1],
0123 (float) global[Acts::eFreePos0], (float) global[Acts::eFreePos1], (float) global[Acts::eFreePos2],
0124 (float) fullCalibrated[Acts::eBoundLoc0], (float) fullCalibrated[Acts::eBoundLoc1]};
0125 tree->Fill(data);
0126
0127 if (chi2 > chi2cut)
0128 {
0129 outlier = true;
0130 }
0131
0132 if (verbosity > 2)
0133 {
0134 std::cout << "Meas vol id and layer " << volume << ", " << layer
0135 << " and chi2cut "
0136 << chi2cut << " so returning outlier : " << outlier
0137 << std::endl;
0138 }
0139
0140 return outlier;
0141 }
0142 };
0143
0144 #endif
