44 :
IShowerTool(pset.get<fhicl::ParameterSet>(
"BaseTools"))
45 , fVerbose(pset.get<int>(
"Verbose"))
46 , fUsePandoraVertex(pset.get<bool>(
"UsePandoraVertex"))
47 , fInitialTrackSpacePointsInputLabel(pset.get<std::string>(
"InitialTrackSpacePointsInputLabel"))
48 , fShowerStartPositionInputLabel(pset.get<std::string>(
"ShowerStartPositionInputLabel"))
49 , fInitialTrackInputLabel(pset.get<std::string>(
"InitialTrackInputLabel"))
50 , fShowerDirectionOutputLabel(pset.get<std::string>(
"ShowerDirectionOutputLabel"))
54 const art::Ptr<recob::PFParticle>& pfparticle,
62 mf::LogError(
"ShowerTrackSpacePointDirection")
63 <<
"Initial track spacepoints not set" << std::endl;
70 mf::LogError(
"ShowerTrackSpacePointDirection")
71 <<
"Start position not set, returning " << std::endl;
76 geo::Point_t StartPosition = {-999, -999, -999};
84 mf::LogError(
"ShowerTrackSpacePointDirection") <<
"Initial track not set" << std::endl;
87 recob::Track InitialTrack;
89 geo::Point_t Start_point = InitialTrack.Start();
90 StartPosition = {Start_point.X(), Start_point.Y(), Start_point.Z()};
94 std::vector<art::Ptr<recob::SpacePoint>> intitaltrack_sp;
98 float sumX = 0, sumX2 = 0;
99 float sumY = 0, sumY2 = 0;
100 float sumZ = 0, sumZ2 = 0;
103 for (
auto const& sp : intitaltrack_sp) {
106 auto const pos = sp->position() - StartPosition;
107 if (pos.R() == 0) {
continue; }
110 sumX2 += pos.X() * pos.X();
112 sumY2 += pos.Y() * pos.Y();
114 sumZ2 += pos.Z() * pos.Z();
117 float NumSps = intitaltrack_sp.size();
118 auto const Mean = geo::Vector_t{sumX / NumSps, sumY / NumSps, sumZ / NumSps}.Unit();
123 if (sumX2 / NumSps - ((sumX / NumSps) * ((sumX / NumSps))) > 0) {
124 RMSX = std::sqrt(sumX2 / NumSps - ((sumX / NumSps) * ((sumX / NumSps))));
126 if (sumY2 / NumSps - ((sumY / NumSps) * ((sumY / NumSps))) > 0) {
127 RMSY = std::sqrt(sumY2 / NumSps - ((sumY / NumSps) * ((sumY / NumSps))));
129 if (sumZ2 / NumSps - ((sumZ / NumSps) * ((sumZ / NumSps))) > 0) {
130 RMSZ = std::sqrt(sumZ2 / NumSps - ((sumZ / NumSps) * ((sumZ / NumSps))));
134 geo::Vector_t Direction_Mean = {0, 0, 0};
136 for (
auto const sp : intitaltrack_sp) {
137 auto const Direction = sp->position() - StartPosition;
138 if ((std::abs((Direction - Mean).X()) < 1 * RMSX) &&
139 (std::abs((Direction - Mean).Y()) < 1 * RMSY) &&
140 (std::abs((Direction - Mean).Z()) < 1 * RMSZ)) {
141 if (Direction.R() == 0) {
continue; }
143 Direction_Mean += Direction;
149 Direction_Mean = Direction_Mean.Unit();
154 mf::LogError(
"ShowerTrackSpacePointDirection")
155 <<
"None of the points are within 1 sigma" << std::endl;