All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 
17 namespace 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
48  ShowerTrackDirection::CalculateElement(const art::Ptr<recob::PFParticle>& pfparticle,
49  art::Event& Event,
50  reco::shower::ShowerElementHolder& ShowerEleHolder)
51  {
52 
53  //Check the Track has been defined
54  if (!ShowerEleHolder.CheckElement("InitialTrack")) {
55  if (fVerbose) mf::LogError("ShowerTrackDirection") << "Initial track not set" << std::endl;
56  return 0;
57  }
58 
59  //Check the start position is set.
60  if (fUsePandoraVertex && !ShowerEleHolder.CheckElement("ShowerStartPosition")) {
61  if (fVerbose)
62  mf::LogError("ShowerTrackDirection") << "Start position not set, returning " << std::endl;
63  return 0;
64  }
65 
66  //Get the track
67  recob::Track InitialTrack;
68  ShowerEleHolder.GetElement("InitialTrack", InitialTrack);
69 
70  if (fUsePositionInfo) {
71  geo::Point_t StartPosition;
72  if (fUsePandoraVertex) {
73  TVector3 StartPosition_vec = {-999, -999, -999};
74  ShowerEleHolder.GetElement("ShowerStartPosition", StartPosition_vec);
75  StartPosition.SetCoordinates(
76  StartPosition_vec.X(), StartPosition_vec.Y(), StartPosition_vec.Z());
77  }
78  else {
79  StartPosition = InitialTrack.Start();
80  }
81 
82  //Calculate the mean direction and the the standard deviation
83  float sumX = 0, sumX2 = 0;
84  float sumY = 0, sumY2 = 0;
85  float sumZ = 0, sumZ2 = 0;
86  for (unsigned int traj = 0; traj < InitialTrack.NumberTrajectoryPoints(); ++traj) {
87 
88  //Ignore bogus flags.
89  auto flags = InitialTrack.FlagsAtPoint(traj);
90  if (flags.isSet(recob::TrajectoryPointFlags::InvalidHitIndex) ||
91  flags.isSet(recob::TrajectoryPointFlagTraits::NoPoint)) {
92  continue;
93  }
94 
95  //Get the direction to the trajectory position.
96  geo::Vector_t TrajPosition = (InitialTrack.LocationAtPoint(traj) - StartPosition).Unit();
97  sumX += TrajPosition.X();
98  sumX2 += TrajPosition.X() * TrajPosition.X();
99  sumY += TrajPosition.Y();
100  sumY2 += TrajPosition.Y() * TrajPosition.Y();
101  sumZ += TrajPosition.Z();
102  sumZ2 += TrajPosition.Z() * TrajPosition.Z();
103  }
104 
105  float NumTraj = InitialTrack.NumberTrajectoryPoints();
106  geo::Vector_t Mean = {sumX / NumTraj, sumY / NumTraj, sumZ / NumTraj};
107  Mean = Mean.Unit();
108 
109  float RMSX = 999;
110  float RMSY = 999;
111  float RMSZ = 999;
112  if (sumX2 / NumTraj - ((sumX / NumTraj) * ((sumX / NumTraj))) > 0) {
113  RMSX = std::sqrt(sumX2 / NumTraj - ((sumX / NumTraj) * ((sumX / NumTraj))));
114  }
115  if (sumY2 / NumTraj - ((sumY / NumTraj) * ((sumY / NumTraj))) > 0) {
116  RMSY = std::sqrt(sumY2 / NumTraj - ((sumY / NumTraj) * ((sumY / NumTraj))));
117  }
118  if (sumZ2 / NumTraj - ((sumZ / NumTraj) * ((sumZ / NumTraj))) > 0) {
119  RMSZ = std::sqrt(sumZ2 / NumTraj - ((sumZ / NumTraj) * ((sumZ / NumTraj))));
120  }
121 
122  TVector3 Direction_Mean = {0, 0, 0};
123  int N = 0;
124  //Remove trajectory points from the mean that are not with one sigma.
125  for (unsigned int traj = 0; traj < InitialTrack.NumberTrajectoryPoints(); ++traj) {
126 
127  //Ignore bogus flags.
128  auto flags = InitialTrack.FlagsAtPoint(traj);
129  if (flags.isSet(recob::TrajectoryPointFlagTraits::NoPoint)) { continue; }
130 
131  //Get the direction of the trajectory point.
132  geo::Point_t TrajPosition = InitialTrack.LocationAtPoint(traj);
133  geo::Vector_t Direction = (TrajPosition - StartPosition).Unit();
134 
135  //Remove points not within 1RMS.
136  if ((std::abs((Direction - Mean).X()) < 1 * RMSX) &&
137  (std::abs((Direction - Mean).Y()) < 1 * RMSY) &&
138  (std::abs((Direction - Mean).Z()) < 1 * RMSZ)) {
139  TVector3 Direction_vec = {Direction.X(), Direction.Y(), Direction.Z()};
140  if (Direction_vec.Mag() == 0) { continue; }
141  Direction_Mean += Direction_vec;
142  ++N;
143  }
144  }
145 
146  //Take the mean value
147  if (N > 0) {
148  TVector3 Direction = Direction_Mean.Unit();
149  TVector3 DirectionErr = {RMSX, RMSY, RMSZ};
150  ShowerEleHolder.SetElement(Direction, DirectionErr, "ShowerDirection");
151  }
152  else {
153  if (fVerbose)
154  mf::LogError("ShowerTrackDirection")
155  << "None of the points are within 1 sigma" << std::endl;
156  return 1;
157  }
158 
159  return 0;
160  }
161  else { // if(fUsePositionInfo)
162 
163  float sumX = 0, sumX2 = 0;
164  float sumY = 0, sumY2 = 0;
165  float sumZ = 0, sumZ2 = 0;
166  for (unsigned int traj = 0; traj < InitialTrack.NumberTrajectoryPoints(); ++traj) {
167 
168  //Ignore bogus points
169  auto flags = InitialTrack.FlagsAtPoint(traj);
170  if (flags.isSet(recob::TrajectoryPointFlagTraits::NoPoint)) { continue; }
171 
172  //Get the direction.
173  geo::Vector_t Direction = InitialTrack.DirectionAtPoint(traj);
174  sumX += Direction.X();
175  sumX2 += Direction.X() * Direction.X();
176  sumY += Direction.Y();
177  sumY2 += Direction.Y() * Direction.Y();
178  sumZ += Direction.Z();
179  sumZ2 += Direction.Z() * Direction.Z();
180  }
181 
182  float NumTraj = InitialTrack.NumberTrajectoryPoints();
183  geo::Vector_t Mean = {sumX / NumTraj, sumY / NumTraj, sumZ / NumTraj};
184  Mean = Mean.Unit();
185 
186  float RMSX = 999;
187  float RMSY = 999;
188  float RMSZ = 999;
189  if (sumX2 / NumTraj - ((sumX / NumTraj) * ((sumX / NumTraj))) > 0) {
190  RMSX = std::sqrt(sumX2 / NumTraj - ((sumX / NumTraj) * ((sumX / NumTraj))));
191  }
192  if (sumY2 / NumTraj - ((sumY / NumTraj) * ((sumY / NumTraj))) > 0) {
193  RMSY = std::sqrt(sumY2 / NumTraj - ((sumY / NumTraj) * ((sumY / NumTraj))));
194  }
195  if (sumZ2 / NumTraj - ((sumZ / NumTraj) * ((sumZ / NumTraj))) > 0) {
196  RMSZ = std::sqrt(sumZ2 / NumTraj - ((sumZ / NumTraj) * ((sumZ / NumTraj))));
197  }
198 
199  //Remove trajectory points from the mean that are not with one sigma.
200  float N = 0.;
201  TVector3 Direction_Mean = {0, 0, 0};
202  for (unsigned int traj = 0; traj < InitialTrack.NumberTrajectoryPoints(); ++traj) {
203 
204  auto flags = InitialTrack.FlagsAtPoint(traj);
205  if (flags.isSet(recob::TrajectoryPointFlagTraits::NoPoint)) { continue; }
206 
207  geo::Vector_t Direction = InitialTrack.DirectionAtPoint(traj).Unit();
208  if ((std::abs((Direction - Mean).X()) < 1 * RMSX) &&
209  (std::abs((Direction - Mean).Y()) < 1 * RMSY) &&
210  (std::abs((Direction - Mean).Z()) < 1 * RMSZ)) {
211  TVector3 Direction_vec = {Direction.X(), Direction.Y(), Direction.Z()};
212  if (Direction_vec.Mag() == 0) { continue; }
213  Direction_Mean += Direction_vec;
214  ++N;
215  }
216  }
217 
218  //Take the mean value
219  if (N > 0) {
220  TVector3 Direction = Direction_Mean.Unit();
221  TVector3 DirectionErr = {RMSX, RMSY, RMSZ};
222  ShowerEleHolder.SetElement(Direction, DirectionErr, "ShowerDirection");
223  }
224  else {
225  if (fVerbose)
226  mf::LogError("ShowerTrackDirection")
227  << "None of the points are within 1 sigma" << std::endl;
228  return 1;
229  }
230  }
231  return 0;
232  }
233 }
234 
235 DEFINE_ART_CLASS_TOOL(ShowerRecoTools::ShowerTrackDirection)
geo::Vector_t
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
recob::TrajectoryPointFlagTraits::NoPoint
static constexpr Flag_t NoPoint
The trajectory point is not defined.
Definition: TrajectoryPointFlags.h:119
recob::Track::LocationAtPoint
Point_t const & LocationAtPoint(size_t i) const
Definition: bj/lardataobj/RecoBase/Track.h:126
get
get()
Definition: PMTconfigurationExtraction_module.cc:168
reco::shower::ShowerElementHolder::SetElement
void SetElement(T &dataproduct, const std::string &Name, bool checktag=false)
Definition: ShowerElementHolder.hh:362
recob::Track::NumberTrajectoryPoints
size_t NumberTrajectoryPoints() const
Various functions related to the presence and the number of (valid) points.
Definition: bj/lardataobj/RecoBase/Track.h:102
larg4::bool
static constexpr bool
Definition: OpFastScintillation.hh:588
X
then echo echo For and will not be changed by echo further linking echo echo B echo The symbol is in the uninitialized data multiple common symbols may appear with the echo same name If the symbol is defined the common echo symbols are treated as undefined references For more echo details on common see the discussion of warn common echo in *Note Linker see the discussion of warn common echo in *Note Linker such as a global int variable echo as opposed to a large global array echo echo I echo The symbol is an indirect reference to another symbol This echo is a GNU extension to the a out object file format which is echo rarely used echo echo N echo The symbol is a debugging symbol echo echo R echo The symbol is in a read only data section echo echo S echo The symbol is in an uninitialized data section for small echo objects echo echo T echo The symbol is in the the normal defined echo symbol is used with no error When a weak undefined symbol echo is linked and the symbol is not the value of the echo weak symbol becomes zero with no error echo echo W echo The symbol is a weak symbol that has not been specifically echo tagged as a weak object symbol When a weak defined symbol echo is linked with a normal defined the normal defined echo symbol is used with no error When a weak undefined symbol echo is linked and the symbol is not the value of the echo weak symbol becomes zero with no error echo echo echo The symbol is a stabs symbol in an a out object file In echo this the next values printed are the stabs other echo the stabs desc and the stab type Stabs symbols are echo used to hold debugging information For more echo see *Note or object file format specific echo echo For Mac OS X
Definition: find_global_symbol.sh:107
ShowerRecoTools::ShowerTrackDirection::ShowerTrackDirection
ShowerTrackDirection(const fhicl::ParameterSet &pset)
Definition: ShowerTrackDirection_tool.cc:40
abs
T abs(T value)
Definition: sparse_vector_test.cc:30
recob::Track::Start
Point_t const & Start() const
Access to track position at different points.
Definition: bj/lardataobj/RecoBase/Track.h:123
reco::shower::ShowerElementHolder::CheckElement
bool CheckElement(const std::string &Name) const
Definition: ShowerElementHolder.hh:417
reco::shower::ShowerElementHolder::GetElement
int GetElement(const std::string &Name, T &Element) const
Definition: ShowerElementHolder.hh:221
recob::TrajectoryPointFlags::InvalidHitIndex
static constexpr HitIndex_t InvalidHitIndex
Value marking an invalid hit index.
Definition: TrajectoryPointFlags.h:336
Track.h
Provides recob::Track data product.
N
process_name largeant stream1 can override from command line with o or output physics producers generator N
Definition: prova_source.fcl:33
recob::Track::FlagsAtPoint
PointFlags_t const & FlagsAtPoint(size_t i) const
Definition: bj/lardataobj/RecoBase/Track.h:118
ShowerRecoTools::ShowerTrackDirection::CalculateElement
int CalculateElement(const art::Ptr< recob::PFParticle > &pfparticle, art::Event &Event, reco::shower::ShowerElementHolder &ShowerEleHolder) override
Definition: ShowerTrackDirection_tool.cc:48
recob::Track::DirectionAtPoint
Vector_t DirectionAtPoint(size_t i) const
Definition: bj/lardataobj/RecoBase/Track.h:134
geo::Point_t
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
recob::Track
Track from a non-cascading particle.A recob::Track consists of a recob::TrackTrajectory, plus additional members relevant for a &quot;fitted&quot; track:
Definition: bj/lardataobj/RecoBase/Track.h:49
Generated by   doxygen 1.8.5