Pandora
Pandora source code navigator
Loading...
Searching...
No Matches
ShowerTrackDirection_tool.cc
Go to the documentation of this file.
1//############################################################################
2//### Name: ShowerTrackDirection ###
3//### Author: Dominic Barker ###
4//### Date: 13.05.19 ###
5//### Description: Tool for finding the shower direction using the ###
6//### initial track method. ###
7//############################################################################
8
9//Framework Includes
10#include "art/Utilities/ToolMacros.h"
11
12//LArSoft Includes
13#include "larpandora/LArPandoraEventBuilding/LArPandoraShower/Tools/IShowerTool.h"
14
15#include "lardataobj/RecoBase/Track.h"
16
17namespace ShowerRecoTools {
18
19 class ShowerTrackDirection : IShowerTool {
20
21 public:
22 ShowerTrackDirection(const fhicl::ParameterSet& pset);
23
24 //Find Track Direction using initial track.
25 int CalculateElement(const art::Ptr<recob::PFParticle>& pfparticle,
26 art::Event& Event,
27 reco::shower::ShowerElementHolder& ShowerEleHolder) override;
28
29 private:
30 //fcl
31 int fVerbose;
32 bool fUsePandoraVertex; //Direction from point defined as
33 //(Position of traj point - Vertex) rather than
34 //(Position of traj point - Track Start Point).
35 bool fUsePositionInfo; //Don't use the directionAtPoint rather
36 //than definition above.
37 //i.e((Position of traj point + 1) - (Position of traj point)
38 };
39
40 ShowerTrackDirection::ShowerTrackDirection(const fhicl::ParameterSet& pset)
41 : IShowerTool(pset.get<fhicl::ParameterSet>("BaseTools"))
42 , fVerbose(pset.get<int>("Verbose"))
43 , fUsePandoraVertex(pset.get<bool>("UsePandoraVertex"))
44 , fUsePositionInfo(pset.get<bool>("UsePositionInfo"))
45 {}
46
47 int ShowerTrackDirection::CalculateElement(const art::Ptr<recob::PFParticle>& pfparticle,
48 art::Event& Event,
49 reco::shower::ShowerElementHolder& ShowerEleHolder)
50 {
51
52 //Check the Track has been defined
53 if (!ShowerEleHolder.CheckElement("InitialTrack")) {
54 if (fVerbose) mf::LogError("ShowerTrackDirection") << "Initial track not set" << std::endl;
55 return 0;
56 }
57
58 //Check the start position is set.
59 if (fUsePandoraVertex && !ShowerEleHolder.CheckElement("ShowerStartPosition")) {
60 if (fVerbose)
61 mf::LogError("ShowerTrackDirection") << "Start position not set, returning " << std::endl;
62 return 0;
63 }
64
65 //Get the track
66 recob::Track InitialTrack;
67 ShowerEleHolder.GetElement("InitialTrack", InitialTrack);
68
69 if (fUsePositionInfo) {
70 geo::Point_t StartPosition;
71 if (fUsePandoraVertex) { ShowerEleHolder.GetElement("ShowerStartPosition", StartPosition); }
72 else {
73 StartPosition = InitialTrack.Start();
74 }
75
76 //Calculate the mean direction and the the standard deviation
77 float sumX = 0, sumX2 = 0;
78 float sumY = 0, sumY2 = 0;
79 float sumZ = 0, sumZ2 = 0;
80 for (unsigned int traj = 0; traj < InitialTrack.NumberTrajectoryPoints(); ++traj) {
81
82 //Ignore bogus flags.
83 auto flags = InitialTrack.FlagsAtPoint(traj);
84 if (flags.isSet(recob::TrajectoryPointFlags::InvalidHitIndex) ||
85 flags.isSet(recob::TrajectoryPointFlagTraits::NoPoint)) {
86 continue;
87 }
88
89 //Get the direction to the trajectory position.
90 geo::Vector_t TrajPosition = (InitialTrack.LocationAtPoint(traj) - StartPosition).Unit();
91 sumX += TrajPosition.X();
92 sumX2 += TrajPosition.X() * TrajPosition.X();
93 sumY += TrajPosition.Y();
94 sumY2 += TrajPosition.Y() * TrajPosition.Y();
95 sumZ += TrajPosition.Z();
96 sumZ2 += TrajPosition.Z() * TrajPosition.Z();
97 }
98
99 float NumTraj = InitialTrack.NumberTrajectoryPoints();
100 geo::Vector_t Mean = {sumX / NumTraj, sumY / NumTraj, sumZ / NumTraj};
101 Mean = Mean.Unit();
102
103 float RMSX = 999;
104 float RMSY = 999;
105 float RMSZ = 999;
106 if (sumX2 / NumTraj - ((sumX / NumTraj) * ((sumX / NumTraj))) > 0) {
107 RMSX = std::sqrt(sumX2 / NumTraj - ((sumX / NumTraj) * ((sumX / NumTraj))));
108 }
109 if (sumY2 / NumTraj - ((sumY / NumTraj) * ((sumY / NumTraj))) > 0) {
110 RMSY = std::sqrt(sumY2 / NumTraj - ((sumY / NumTraj) * ((sumY / NumTraj))));
111 }
112 if (sumZ2 / NumTraj - ((sumZ / NumTraj) * ((sumZ / NumTraj))) > 0) {
113 RMSZ = std::sqrt(sumZ2 / NumTraj - ((sumZ / NumTraj) * ((sumZ / NumTraj))));
114 }
115
116 TVector3 Direction_Mean = {0, 0, 0};
117 int N = 0;
118 //Remove trajectory points from the mean that are not with one sigma.
119 for (unsigned int traj = 0; traj < InitialTrack.NumberTrajectoryPoints(); ++traj) {
120
121 //Ignore bogus flags.
122 auto flags = InitialTrack.FlagsAtPoint(traj);
123 if (flags.isSet(recob::TrajectoryPointFlagTraits::NoPoint)) { continue; }
124
125 //Get the direction of the trajectory point.
126 geo::Point_t TrajPosition = InitialTrack.LocationAtPoint(traj);
127 geo::Vector_t Direction = (TrajPosition - StartPosition).Unit();
128
129 //Remove points not within 1RMS.
130 if (auto MeanSubtractedDir = Direction - Mean; (std::abs(MeanSubtractedDir.X()) < RMSX) &&
131 (std::abs(MeanSubtractedDir.Y()) < RMSY) &&
132 (std::abs(MeanSubtractedDir.Z()) < RMSZ)) {
133 if (Direction.R() == 0) { continue; }
134 Direction_Mean += geo::vect::convertTo<TVector3>(Direction);
135 ++N;
136 }
137 }
138
139 //Take the mean value
140 if (N > 0) {
141 geo::Vector_t Direction = geo::vect::toVector(Direction_Mean.Unit());
142 geo::Vector_t DirectionErr = {RMSX, RMSY, RMSZ};
143 ShowerEleHolder.SetElement(Direction, DirectionErr, "ShowerDirection");
144 }
145 else {
146 if (fVerbose)
147 mf::LogError("ShowerTrackDirection")
148 << "None of the points are within 1 sigma" << std::endl;
149 return 1;
150 }
151
152 return 0;
153 }
154 else { // if(fUsePositionInfo)
155
156 float sumX = 0, sumX2 = 0;
157 float sumY = 0, sumY2 = 0;
158 float sumZ = 0, sumZ2 = 0;
159 for (unsigned int traj = 0; traj < InitialTrack.NumberTrajectoryPoints(); ++traj) {
160
161 //Ignore bogus points
162 auto flags = InitialTrack.FlagsAtPoint(traj);
163 if (flags.isSet(recob::TrajectoryPointFlagTraits::NoPoint)) { continue; }
164
165 //Get the direction.
166 geo::Vector_t Direction = InitialTrack.DirectionAtPoint(traj);
167 sumX += Direction.X();
168 sumX2 += Direction.X() * Direction.X();
169 sumY += Direction.Y();
170 sumY2 += Direction.Y() * Direction.Y();
171 sumZ += Direction.Z();
172 sumZ2 += Direction.Z() * Direction.Z();
173 }
174
175 float NumTraj = InitialTrack.NumberTrajectoryPoints();
176 geo::Vector_t Mean = {sumX / NumTraj, sumY / NumTraj, sumZ / NumTraj};
177 Mean = Mean.Unit();
178
179 float RMSX = 999;
180 float RMSY = 999;
181 float RMSZ = 999;
182 if (sumX2 / NumTraj - ((sumX / NumTraj) * ((sumX / NumTraj))) > 0) {
183 RMSX = std::sqrt(sumX2 / NumTraj - ((sumX / NumTraj) * ((sumX / NumTraj))));
184 }
185 if (sumY2 / NumTraj - ((sumY / NumTraj) * ((sumY / NumTraj))) > 0) {
186 RMSY = std::sqrt(sumY2 / NumTraj - ((sumY / NumTraj) * ((sumY / NumTraj))));
187 }
188 if (sumZ2 / NumTraj - ((sumZ / NumTraj) * ((sumZ / NumTraj))) > 0) {
189 RMSZ = std::sqrt(sumZ2 / NumTraj - ((sumZ / NumTraj) * ((sumZ / NumTraj))));
190 }
191
192 //Remove trajectory points from the mean that are not with one sigma.
193 float N = 0.;
194 TVector3 Direction_Mean = {0, 0, 0};
195 for (unsigned int traj = 0; traj < InitialTrack.NumberTrajectoryPoints(); ++traj) {
196
197 auto flags = InitialTrack.FlagsAtPoint(traj);
198 if (flags.isSet(recob::TrajectoryPointFlagTraits::NoPoint)) { continue; }
199
200 geo::Vector_t Direction = InitialTrack.DirectionAtPoint(traj).Unit();
201 if ((std::abs((Direction - Mean).X()) < 1 * RMSX) &&
202 (std::abs((Direction - Mean).Y()) < 1 * RMSY) &&
203 (std::abs((Direction - Mean).Z()) < 1 * RMSZ)) {
204 TVector3 Direction_vec = {Direction.X(), Direction.Y(), Direction.Z()};
205 if (Direction_vec.Mag() == 0) { continue; }
206 Direction_Mean += Direction_vec;
207 ++N;
208 }
209 }
210
211 //Take the mean value
212 if (N > 0) {
213 geo::Vector_t Direction = geo::vect::toVector(Direction_Mean.Unit());
214 geo::Vector_t DirectionErr = {RMSX, RMSY, RMSZ};
215 ShowerEleHolder.SetElement(Direction, DirectionErr, "ShowerDirection");
216 }
217 else {
218 if (fVerbose)
219 mf::LogError("ShowerTrackDirection")
220 << "None of the points are within 1 sigma" << std::endl;
221 return 1;
222 }
223 }
224 return 0;
225 }
226}
227
228DEFINE_ART_CLASS_TOOL(ShowerRecoTools::ShowerTrackDirection)
ShowerRecoTools::ShowerTrackDirection::ShowerTrackDirection
ShowerTrackDirection(const fhicl::ParameterSet &pset)
Definition ShowerTrackDirection_tool.cc:40
ShowerRecoTools::ShowerTrackDirection::CalculateElement
int CalculateElement(const art::Ptr< recob::PFParticle > &pfparticle, art::Event &Event, reco::shower::ShowerElementHolder &ShowerEleHolder) override
Definition ShowerTrackDirection_tool.cc:47
reco::shower::ShowerElementHolder::GetElement
int GetElement(const std::string &Name, T &Element) const
Definition ShowerElementHolder.hh:219
reco::shower::ShowerElementHolder::SetElement
void SetElement(T &dataproduct, const std::string &Name, bool checktag=false)
Definition ShowerElementHolder.hh:405
reco::shower::ShowerElementHolder::CheckElement
bool CheckElement(const std::string &Name) const
Definition ShowerElementHolder.hh:470