Pandora
Pandora source code navigator
Loading...
Searching...
No Matches
ShowerTrackSpacePointDirection_tool.cc
Go to the documentation of this file.
1//############################################################################
2//### Name: ShowerTrackSpacePointDirection ###
3//### Author: Dominic Barker ###
4//### Date: 13.05.19 ###
5//### Description: Tool for finding the shower direction using the ###
6//### the average direction of theinitial track spacepoints. ###
7//############################################################################
8
9//Framework Includes
10#include "art/Utilities/ToolMacros.h"
11
12//LArSoft Includes
14
15#include "lardataobj/RecoBase/SpacePoint.h"
16#include "lardataobj/RecoBase/Track.h"
17
18namespace ShowerRecoTools {
19
21
22 public:
23 ShowerTrackSpacePointDirection(const fhicl::ParameterSet& pset);
24
25 //Calculate the direction using the initial track spacepoints
26 int CalculateElement(const art::Ptr<recob::PFParticle>& pfparticle,
27 art::Event& Event,
28 reco::shower::ShowerElementHolder& ShowerEleHolder) override;
29
30 private:
31 //fcl
33 bool fUsePandoraVertex; //Direction from point defined as
34 //(Position of SP - Vertex) rather than
35 //(Position of SP - Track Start Point).
36
41 };
42
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"))
51 {}
52
54 const art::Ptr<recob::PFParticle>& pfparticle,
55 art::Event& Event,
56 reco::shower::ShowerElementHolder& ShowerEleHolder)
57 {
58
59 //Check the Track Hits has been defined
61 if (fVerbose)
62 mf::LogError("ShowerTrackSpacePointDirection")
63 << "Initial track spacepoints not set" << std::endl;
64 return 0;
65 }
66
67 //Check the start position is set.
69 if (fVerbose)
70 mf::LogError("ShowerTrackSpacePointDirection")
71 << "Start position not set, returning " << std::endl;
72 return 0;
73 }
74
75 //Get the start poistion
76 geo::Point_t StartPosition = {-999, -999, -999};
78 ShowerEleHolder.GetElement(fShowerStartPositionInputLabel, StartPosition);
79 }
80 else {
81 //Check the Tracks has been defined
82 if (!ShowerEleHolder.CheckElement(fInitialTrackInputLabel)) {
83 if (fVerbose)
84 mf::LogError("ShowerTrackSpacePointDirection") << "Initial track not set" << std::endl;
85 return 0;
86 }
87 recob::Track InitialTrack;
88 ShowerEleHolder.GetElement(fInitialTrackInputLabel, InitialTrack);
89 geo::Point_t Start_point = InitialTrack.Start();
90 StartPosition = {Start_point.X(), Start_point.Y(), Start_point.Z()};
91 }
92
93 //Get the initial track hits.
94 std::vector<art::Ptr<recob::SpacePoint>> intitaltrack_sp;
95 ShowerEleHolder.GetElement(fInitialTrackSpacePointsInputLabel, intitaltrack_sp);
96
97 //Calculate the mean direction and the the standard deviation
98 float sumX = 0, sumX2 = 0;
99 float sumY = 0, sumY2 = 0;
100 float sumZ = 0, sumZ2 = 0;
101
102 //Get the spacepoints associated to the track hit
103 for (auto const& sp : intitaltrack_sp) {
104
105 //Get the direction relative to the start positon
106 auto const pos = sp->position() - StartPosition;
107 if (pos.R() == 0) { continue; }
108
109 sumX = pos.X();
110 sumX2 += pos.X() * pos.X();
111 sumY = pos.Y();
112 sumY2 += pos.Y() * pos.Y();
113 sumZ = pos.Z();
114 sumZ2 += pos.Z() * pos.Z();
115 }
116
117 float NumSps = intitaltrack_sp.size();
118 auto const Mean = geo::Vector_t{sumX / NumSps, sumY / NumSps, sumZ / NumSps}.Unit();
119
120 float RMSX = 999;
121 float RMSY = 999;
122 float RMSZ = 999;
123 if (sumX2 / NumSps - ((sumX / NumSps) * ((sumX / NumSps))) > 0) {
124 RMSX = std::sqrt(sumX2 / NumSps - ((sumX / NumSps) * ((sumX / NumSps))));
125 }
126 if (sumY2 / NumSps - ((sumY / NumSps) * ((sumY / NumSps))) > 0) {
127 RMSY = std::sqrt(sumY2 / NumSps - ((sumY / NumSps) * ((sumY / NumSps))));
128 }
129 if (sumZ2 / NumSps - ((sumZ / NumSps) * ((sumZ / NumSps))) > 0) {
130 RMSZ = std::sqrt(sumZ2 / NumSps - ((sumZ / NumSps) * ((sumZ / NumSps))));
131 }
132
133 //Loop over the spacepoints and remove ones the relative direction is not within one sigma.
134 geo::Vector_t Direction_Mean = {0, 0, 0};
135 int N = 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; }
142 ++N;
143 Direction_Mean += Direction;
144 }
145 }
146
147 if (N > 0) {
148 //Take the mean value
149 Direction_Mean = Direction_Mean.Unit();
150 ShowerEleHolder.SetElement(Direction_Mean, fShowerDirectionOutputLabel);
151 }
152 else {
153 if (fVerbose)
154 mf::LogError("ShowerTrackSpacePointDirection")
155 << "None of the points are within 1 sigma" << std::endl;
156 return 1;
157 }
158 return 0;
159 }
160}
161
int CalculateElement(const art::Ptr< recob::PFParticle > &pfparticle, art::Event &Event, reco::shower::ShowerElementHolder &ShowerEleHolder) override
int GetElement(const std::string &Name, T &Element) const
void SetElement(T &dataproduct, const std::string &Name, bool checktag=false)
bool CheckElement(const std::string &Name) const