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
11 
12 //LArSoft Includes
14 
15 namespace ShowerRecoTools {
16 
18 
19  public:
21 
22  //Find Track Direction using initial track.
23  int CalculateElement(const art::Ptr<recob::PFParticle>& pfparticle,
25  reco::shower::ShowerElementHolder& ShowerEleHolder) override;
26 
27  private:
28  //fcl
29  int fVerbose;
30  bool fUsePandoraVertex; //Direction from point defined as
31  //(Position of traj point - Vertex) rather than
32  //(Position of traj point - Track Start Point).
33  bool fUsePositionInfo; //Don't use the directionAtPoint rather
34  //than definition above.
35  //i.e((Position of traj point + 1) - (Position of traj point)
36  };
37 
39  : IShowerTool(pset.get<fhicl::ParameterSet>("BaseTools"))
40  , fVerbose(pset.get<int>("Verbose"))
41  , fUsePandoraVertex(pset.get<bool>("UsePandoraVertex"))
42  , fUsePositionInfo(pset.get<bool>("UsePositionInfo"))
43  {}
44 
45  int
47  art::Event& Event,
48  reco::shower::ShowerElementHolder& ShowerEleHolder)
49  {
50 
51  //Check the Track has been defined
52  if (!ShowerEleHolder.CheckElement("InitialTrack")) {
53  if (fVerbose) mf::LogError("ShowerTrackDirection") << "Initial track not set" << std::endl;
54  return 0;
55  }
56 
57  //Check the start position is set.
58  if (fUsePandoraVertex && !ShowerEleHolder.CheckElement("ShowerStartPosition")) {
59  if (fVerbose)
60  mf::LogError("ShowerTrackDirection") << "Start position not set, returning " << std::endl;
61  return 0;
62  }
63 
64  //Get the track
65  recob::Track InitialTrack;
66  ShowerEleHolder.GetElement("InitialTrack", InitialTrack);
67 
68  if (fUsePositionInfo) {
69  geo::Point_t StartPosition;
70  if (fUsePandoraVertex) {
71  TVector3 StartPosition_vec = {-999, -999, -999};
72  ShowerEleHolder.GetElement("ShowerStartPosition", StartPosition_vec);
73  StartPosition.SetCoordinates(
74  StartPosition_vec.X(), StartPosition_vec.Y(), StartPosition_vec.Z());
75  }
76  else {
77  StartPosition = InitialTrack.Start();
78  }
79 
80  //Calculate the mean direction and the the standard deviation
81  float sumX = 0, sumX2 = 0;
82  float sumY = 0, sumY2 = 0;
83  float sumZ = 0, sumZ2 = 0;
84  for (unsigned int traj = 0; traj < InitialTrack.NumberTrajectoryPoints(); ++traj) {
85 
86  //Ignore bogus flags.
87  auto flags = InitialTrack.FlagsAtPoint(traj);
90  continue;
91  }
92 
93  //Get the direction to the trajectory position.
94  geo::Vector_t TrajPosition = (InitialTrack.LocationAtPoint(traj) - StartPosition).Unit();
95  sumX += TrajPosition.X();
96  sumX2 += TrajPosition.X() * TrajPosition.X();
97  sumY += TrajPosition.Y();
98  sumY2 += TrajPosition.Y() * TrajPosition.Y();
99  sumZ += TrajPosition.Z();
100  sumZ2 += TrajPosition.Z() * TrajPosition.Z();
101  }
102 
103  float NumTraj = InitialTrack.NumberTrajectoryPoints();
104  geo::Vector_t Mean = {sumX / NumTraj, sumY / NumTraj, sumZ / NumTraj};
105  Mean = Mean.Unit();
106 
107  float RMSX = 999;
108  float RMSY = 999;
109  float RMSZ = 999;
110  if (sumX2 / NumTraj - ((sumX / NumTraj) * ((sumX / NumTraj))) > 0) {
111  RMSX = std::sqrt(sumX2 / NumTraj - ((sumX / NumTraj) * ((sumX / NumTraj))));
112  }
113  if (sumY2 / NumTraj - ((sumY / NumTraj) * ((sumY / NumTraj))) > 0) {
114  RMSY = std::sqrt(sumY2 / NumTraj - ((sumY / NumTraj) * ((sumY / NumTraj))));
115  }
116  if (sumZ2 / NumTraj - ((sumZ / NumTraj) * ((sumZ / NumTraj))) > 0) {
117  RMSZ = std::sqrt(sumZ2 / NumTraj - ((sumZ / NumTraj) * ((sumZ / NumTraj))));
118  }
119 
120  TVector3 Direction_Mean = {0, 0, 0};
121  int N = 0;
122  //Remove trajectory points from the mean that are not with one sigma.
123  for (unsigned int traj = 0; traj < InitialTrack.NumberTrajectoryPoints(); ++traj) {
124 
125  //Ignore bogus flags.
126  auto flags = InitialTrack.FlagsAtPoint(traj);
127  if (flags.isSet(recob::TrajectoryPointFlagTraits::NoPoint)) { continue; }
128 
129  //Get the direction of the trajectory point.
130  geo::Point_t TrajPosition = InitialTrack.LocationAtPoint(traj);
131  geo::Vector_t Direction = (TrajPosition - StartPosition).Unit();
132 
133  //Remove points not within 1RMS.
134  if ((std::abs((Direction - Mean).X()) < 1 * RMSX) &&
135  (std::abs((Direction - Mean).Y()) < 1 * RMSY) &&
136  (std::abs((Direction - Mean).Z()) < 1 * RMSZ)) {
137  TVector3 Direction_vec = {Direction.X(), Direction.Y(), Direction.Z()};
138  if (Direction_vec.Mag() == 0) { continue; }
139  Direction_Mean += Direction_vec;
140  ++N;
141  }
142  }
143 
144  //Take the mean value
145  if (N > 0) {
146  TVector3 Direction = Direction_Mean.Unit();
147  TVector3 DirectionErr = {RMSX, RMSY, RMSZ};
148  ShowerEleHolder.SetElement(Direction, DirectionErr, "ShowerDirection");
149  }
150  else {
151  if (fVerbose)
152  mf::LogError("ShowerTrackDirection")
153  << "None of the points are within 1 sigma" << std::endl;
154  return 1;
155  }
156 
157  return 0;
158  }
159  else { // if(fUsePositionInfo)
160 
161  float sumX = 0, sumX2 = 0;
162  float sumY = 0, sumY2 = 0;
163  float sumZ = 0, sumZ2 = 0;
164  for (unsigned int traj = 0; traj < InitialTrack.NumberTrajectoryPoints(); ++traj) {
165 
166  //Ignore bogus points
167  auto flags = InitialTrack.FlagsAtPoint(traj);
168  if (flags.isSet(recob::TrajectoryPointFlagTraits::NoPoint)) { continue; }
169 
170  //Get the direction.
171  geo::Vector_t Direction = InitialTrack.DirectionAtPoint(traj);
172  sumX += Direction.X();
173  sumX2 += Direction.X() * Direction.X();
174  sumY += Direction.Y();
175  sumY2 += Direction.Y() * Direction.Y();
176  sumZ += Direction.Z();
177  sumZ2 += Direction.Z() * Direction.Z();
178  }
179 
180  float NumTraj = InitialTrack.NumberTrajectoryPoints();
181  geo::Vector_t Mean = {sumX / NumTraj, sumY / NumTraj, sumZ / NumTraj};
182  Mean = Mean.Unit();
183 
184  float RMSX = 999;
185  float RMSY = 999;
186  float RMSZ = 999;
187  if (sumX2 / NumTraj - ((sumX / NumTraj) * ((sumX / NumTraj))) > 0) {
188  RMSX = std::sqrt(sumX2 / NumTraj - ((sumX / NumTraj) * ((sumX / NumTraj))));
189  }
190  if (sumY2 / NumTraj - ((sumY / NumTraj) * ((sumY / NumTraj))) > 0) {
191  RMSY = std::sqrt(sumY2 / NumTraj - ((sumY / NumTraj) * ((sumY / NumTraj))));
192  }
193  if (sumZ2 / NumTraj - ((sumZ / NumTraj) * ((sumZ / NumTraj))) > 0) {
194  RMSZ = std::sqrt(sumZ2 / NumTraj - ((sumZ / NumTraj) * ((sumZ / NumTraj))));
195  }
196 
197  //Remove trajectory points from the mean that are not with one sigma.
198  float N = 0.;
199  TVector3 Direction_Mean = {0, 0, 0};
200  for (unsigned int traj = 0; traj < InitialTrack.NumberTrajectoryPoints(); ++traj) {
201 
202  auto flags = InitialTrack.FlagsAtPoint(traj);
203  if (flags.isSet(recob::TrajectoryPointFlagTraits::NoPoint)) { continue; }
204 
205  geo::Vector_t Direction = InitialTrack.DirectionAtPoint(traj).Unit();
206  if ((std::abs((Direction - Mean).X()) < 1 * RMSX) &&
207  (std::abs((Direction - Mean).Y()) < 1 * RMSY) &&
208  (std::abs((Direction - Mean).Z()) < 1 * RMSZ)) {
209  TVector3 Direction_vec = {Direction.X(), Direction.Y(), Direction.Z()};
210  if (Direction_vec.Mag() == 0) { continue; }
211  Direction_Mean += Direction_vec;
212  ++N;
213  }
214  }
215 
216  //Take the mean value
217  if (N > 0) {
218  TVector3 Direction = Direction_Mean.Unit();
219  TVector3 DirectionErr = {RMSX, RMSY, RMSZ};
220  ShowerEleHolder.SetElement(Direction, DirectionErr, "ShowerDirection");
221  }
222  else {
223  if (fVerbose)
224  mf::LogError("ShowerTrackDirection")
225  << "None of the points are within 1 sigma" << std::endl;
226  return 1;
227  }
228  }
229  return 0;
230  }
231 }
232 
#define DEFINE_ART_CLASS_TOOL(tool)
Definition: ToolMacros.h:42
static constexpr Flag_t NoPoint
The trajectory point is not defined.
Point_t const & LocationAtPoint(size_t i) const
Definition: Track.h:126
void SetElement(T &dataproduct, const std::string &Name, bool checktag=false)
size_t NumberTrajectoryPoints() const
Various functions related to the presence and the number of (valid) points.
Definition: Track.h:102
ROOT::Math::DisplacementVector3D< ROOT::Math::Cartesian3D< double >, ROOT::Math::GlobalCoordinateSystemTag > Vector_t
Type for representation of momenta in 3D space.
Definition: geo_vectors.h:164
MaybeLogger_< ELseverityLevel::ELsev_error, false > LogError
ShowerTrackDirection(const fhicl::ParameterSet &pset)
T abs(T value)
Point_t const & Start() const
Access to track position at different points.
Definition: Track.h:123
bool CheckElement(const std::string &Name) const
int GetElement(const std::string &Name, T &Element) const
static constexpr HitIndex_t InvalidHitIndex
Value marking an invalid hit index.
ROOT::Math::PositionVector3D< ROOT::Math::Cartesian3D< double >, ROOT::Math::GlobalCoordinateSystemTag > Point_t
Type for representation of position in physical 3D space.
Definition: geo_vectors.h:184
Definition: types.h:32
PointFlags_t const & FlagsAtPoint(size_t i) const
Definition: Track.h:118
int CalculateElement(const art::Ptr< recob::PFParticle > &pfparticle, art::Event &Event, reco::shower::ShowerElementHolder &ShowerEleHolder) override
Vector_t DirectionAtPoint(size_t i) const
Definition: Track.h:134
Direction
Definition: AssnsIter.h:13
auto const & get(AssnsNode< L, R, D > const &r)
Definition: AssnsNode.h:115
int bool
Definition: qglobal.h:345
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
QTextStream & endl(QTextStream &s)