LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
LArPandoraShowerCheatingAlg.cxx
Go to the documentation of this file.
3 
5 
7 
8 #include "TCanvas.h"
9 #include "TH3.h"
10 #include "TPolyLine3D.h"
11 #include "TPolyMarker3D.h"
12 #include "TString.h"
13 #include "TStyle.h"
14 
15 #include <memory>
16 
18  : fLArPandoraShowerAlg(pset.get<fhicl::ParameterSet>("LArPandoraShowerAlg"))
19  , fHitModuleLabel(pset.get<art::InputTag>("HitModuleLabel"))
20  , fPFParticleLabel(pset.get<art::InputTag>("PFParticleLabel"))
21  , fShowerStartPositionInputLabel(pset.get<std::string>("ShowerStartPositionInputLabel"))
22  , fShowerDirectionInputLabel(pset.get<std::string>("ShowerDirectionInputLabel"))
23  , fInitialTrackSpacePointsInputLabel(pset.get<std::string>("InitialTrackSpacePointsInputLabel"))
24 {}
25 
26 std::map<int, const simb::MCParticle*> shower::LArPandoraShowerCheatingAlg::GetTrueParticleMap()
27  const
28 {
29 
30  const sim::ParticleList& particles = particleInventory->ParticleList();
31 
32  std::map<int, const simb::MCParticle*> trueParticles;
33  // Make a map of track id to particle
34  for (sim::ParticleList::const_iterator particleIt = particles.begin();
35  particleIt != particles.end();
36  ++particleIt) {
37  const simb::MCParticle* particle = particleIt->second;
38  trueParticles[particle->TrackId()] = particle;
39  //std::cout<<"Particle ID: "<<particle->TrackId()<<" and PDG: "<<particle->PdgCode()<<std::endl;
40  }
41  return trueParticles;
42 }
43 
44 std::map<int, std::vector<int>> shower::LArPandoraShowerCheatingAlg::GetTrueChain(
45  std::map<int, const simb::MCParticle*>& trueParticles) const
46 {
47 
48  // Roll up showers if not already done:
49  std::map<int, std::vector<int>> showerMothers;
50 
51  // Loop over daughters and find th`e mothers
52  for (const auto& particleIt : trueParticles) {
53  const simb::MCParticle* particle = particleIt.second;
54  const simb::MCParticle* mother = particle;
55 
56  if (std::abs(particle->PdgCode()) != 11 && std::abs(particle->PdgCode()) != 22) { continue; }
57 
58  // While the grand mother exists and is an electron or photon
59  // Note the true mother will skip this loop and fill itself into the map
60  while (mother->Mother() != 0 && trueParticles.find(mother->Mother()) != trueParticles.end()) {
61 
62  int motherId = mother->Mother();
63  if (std::abs(trueParticles[motherId]->PdgCode()) != 11 &&
64  std::abs(trueParticles[motherId]->PdgCode()) != 22) {
65  break;
66  }
67  mother = trueParticles[motherId];
68  }
69  showerMothers[mother->TrackId()].push_back(particle->TrackId());
70  }
71  return showerMothers;
72 }
73 
75  detinfo::DetectorClocksData const& clockData,
76  const simb::MCParticle* trueParticle,
77  art::Event const& Event,
78  reco::shower::ShowerElementHolder& ShowerEleHolder,
79  const art::Ptr<recob::PFParticle>& pfparticle) const
80 {
81 
82  std::cout << "Making Debug Event Display" << std::endl;
83 
84  //Function for drawing reco showers to check direction and initial track selection
85 
86  // Get run info to make unique canvas names
87  int run = Event.run();
88  int subRun = Event.subRun();
89  int event = Event.event();
90  int PFPID = pfparticle->Self();
91 
92  // Create the canvas
93  TString canvasName = Form("canvas_%i_%i_%i_%i", run, subRun, event, PFPID);
94  TCanvas* canvas = tfs->make<TCanvas>(canvasName, canvasName);
95 
96  // Initialise variables
97  double x = 0;
98  double y = 0;
99  double z = 0;
100 
101  std::vector<art::Ptr<recob::SpacePoint>> showerSpacePoints;
102  std::vector<art::Ptr<recob::SpacePoint>> otherSpacePoints;
103 
104  auto const hitHandle = Event.getValidHandle<std::vector<recob::Hit>>(fHitModuleLabel);
105  std::vector<art::Ptr<recob::Hit>> hits;
106  art::fill_ptr_vector(hits, hitHandle);
107 
108  // Get the hits associated with the space points
109  const art::FindManyP<recob::SpacePoint> fmsph(hitHandle, Event, fPFParticleLabel);
110  if (!fmsph.isValid()) {
111  throw cet::exception("LArPandoraShowerCheatingAlg")
112  << "Spacepoint and hit association not valid. Stopping.";
113  }
114 
115  std::map<art::Ptr<recob::SpacePoint>, art::Ptr<recob::Hit>> spacePointHitMap;
116  //Get the hits from the true particle
117  for (auto hit : hits) {
118  int trueParticleID = std::abs(TrueParticleID(clockData, hit));
119  std::vector<art::Ptr<recob::SpacePoint>> sps = fmsph.at(hit.key());
120  if (sps.size() == 1) {
121  art::Ptr<recob::SpacePoint> sp = sps.front();
122  if (trueParticleID == trueParticle->TrackId()) { showerSpacePoints.push_back(sp); }
123  else {
124  otherSpacePoints.push_back(sp);
125  }
126  }
127  }
128 
129  if (!ShowerEleHolder.CheckElement(fShowerStartPositionInputLabel)) {
130  mf::LogError("LArPandoraShowerCheatingAlg")
131  << "Start position not set, returning " << std::endl;
132  return;
133  }
134  if (!ShowerEleHolder.CheckElement(fShowerDirectionInputLabel)) {
135  mf::LogError("LArPandoraShowerCheatingAlg") << "Direction not set, returning " << std::endl;
136  return;
137  }
138  if (!ShowerEleHolder.CheckElement(fInitialTrackSpacePointsInputLabel)) {
139  mf::LogError("LArPandoraShowerCheatingAlg")
140  << "TrackSpacePoints not set, returning " << std::endl;
141  return;
142  }
143 
144  // Get info from shower property holder
145  geo::Point_t showerStartPosition = {-999, -999, -999};
146  geo::Vector_t showerDirection = {-999, -999, -999};
147  std::vector<art::Ptr<recob::SpacePoint>> trackSpacePoints;
148 
149  ShowerEleHolder.GetElement(fShowerStartPositionInputLabel, showerStartPosition);
150  ShowerEleHolder.GetElement(fShowerDirectionInputLabel, showerDirection);
151  ShowerEleHolder.GetElement(fInitialTrackSpacePointsInputLabel, trackSpacePoints);
152 
153  // Create 3D point at vertex, chosed to be origin for ease of use of display
154  double startXYZ[3] = {0, 0, 0};
155  auto startPoly = std::make_unique<TPolyMarker3D>(1, startXYZ);
156 
157  // Get the min and max projections along the direction to know how long to draw
158  // the direction line
159  double minProj = std::numeric_limits<double>::max();
160  double maxProj = -std::numeric_limits<double>::max();
161 
162  double x_min = std::numeric_limits<double>::max(), x_max = -std::numeric_limits<double>::max();
163  double y_min = std::numeric_limits<double>::max(), y_max = -std::numeric_limits<double>::max();
164  double z_min = std::numeric_limits<double>::max(), z_max = -std::numeric_limits<double>::max();
165 
166  //initialise counter point
167  int point = 0;
168 
169  // Make 3D points for each spacepoint in the shower
170  auto showerPoly = std::make_unique<TPolyMarker3D>(showerSpacePoints.size());
171  for (auto const& spacePoint : showerSpacePoints) {
172  auto const pos = spacePoint->position() - showerStartPosition;
173  x = pos.X();
174  y = pos.Y();
175  z = pos.Z();
176  x_min = std::min(x, x_min);
177  x_max = std::max(x, x_max);
178  y_min = std::min(y, y_min);
179  y_max = std::max(y, y_max);
180  z_min = std::min(z, z_min);
181  z_max = std::max(z, z_max);
182 
183  showerPoly->SetPoint(point, x, y, z);
184  ++point;
185 
186  // Calculate the projection of (point-startpoint) along the direction
187  double proj =
188  fLArPandoraShowerAlg.SpacePointProjection(spacePoint, showerStartPosition, showerDirection);
189 
190  maxProj = std::max(proj, maxProj);
191  minProj = std::min(proj, minProj);
192  } // loop over spacepoints
193 
194  // Create TPolyLine3D arrays
195  double xDirPoints[2] = {minProj * showerDirection.X(), maxProj * showerDirection.X()};
196  double yDirPoints[2] = {minProj * showerDirection.Y(), maxProj * showerDirection.Y()};
197  double zDirPoints[2] = {minProj * showerDirection.Z(), maxProj * showerDirection.Z()};
198 
199  auto dirPoly = std::make_unique<TPolyLine3D>(2, xDirPoints, yDirPoints, zDirPoints);
200 
201  point = 0; // re-initialise counter
202  auto trackPoly = std::make_unique<TPolyMarker3D>(trackSpacePoints.size());
203  for (auto const& spacePoint : trackSpacePoints) {
204  auto const pos = spacePoint->position() - showerStartPosition;
205  x = pos.X();
206  y = pos.Y();
207  z = pos.Z();
208  trackPoly->SetPoint(point, x, y, z);
209  ++point;
210  } // loop over track spacepoints
211 
212  // we want to draw all of the PFParticles in the event
213  //Get the PFParticles
214 
215  auto otherPoly = std::make_unique<TPolyMarker3D>(otherSpacePoints.size());
216 
217  // initialise counters
218  point = 0; // re-initialise counter
219 
220  for (auto const& sp : otherSpacePoints) {
221  auto const pos = sp->position() - showerStartPosition;
222  x = pos.X();
223  y = pos.Y();
224  z = pos.Z();
225  x_min = std::min(x, x_min);
226  x_max = std::max(x, x_max);
227  y_min = std::min(y, y_min);
228  y_max = std::max(y, y_max);
229  z_min = std::min(z, z_min);
230  z_max = std::max(z, z_max);
231  otherPoly->SetPoint(point, x, y, z);
232  ++point;
233  }
234 
235  gStyle->SetOptStat(0);
236  TH3F axes("axes", "", 1, x_min, x_max, 1, y_min, y_max, 1, z_min, z_max);
237  axes.SetDirectory(0);
238  axes.GetXaxis()->SetTitle("X");
239  axes.GetYaxis()->SetTitle("Y");
240  axes.GetZaxis()->SetTitle("Z");
241  axes.Draw();
242 
243  otherPoly->SetMarkerStyle(20);
244  otherPoly->SetMarkerColor(4);
245  otherPoly->Draw();
246 
247  // Draw all of the things
248  showerPoly->SetMarkerStyle(20);
249  showerPoly->Draw();
250  trackPoly->SetMarkerStyle(20);
251  trackPoly->SetMarkerColor(2);
252  trackPoly->Draw();
253  startPoly->SetMarkerStyle(21);
254  startPoly->SetMarkerSize(2);
255  startPoly->SetMarkerColor(3);
256  startPoly->Draw();
257  dirPoly->SetLineWidth(1);
258  dirPoly->SetLineColor(6);
259  dirPoly->Draw();
260 
261  // Save the canvas. Don't usually need this when using TFileService but this in the alg
262  // not a module and didn't work without this so im going with it.
263  canvas->Write();
264 }
265 
267  detinfo::DetectorClocksData const& clockData,
268  const art::Ptr<recob::Hit>& hit) const
269 {
270 
271  double particleEnergy = 0;
272  int likelyTrackID = 0;
274  std::vector<sim::TrackIDE> trackIDs = bt_serv->HitToTrackIDEs(clockData, hit);
275  for (unsigned int idIt = 0; idIt < trackIDs.size(); ++idIt) {
276  if (trackIDs.at(idIt).energy > particleEnergy) {
277  particleEnergy = trackIDs.at(idIt).energy;
278  likelyTrackID = trackIDs.at(idIt).trackID;
279  }
280  }
281  return likelyTrackID;
282 }
283 
285  detinfo::DetectorClocksData const& clockData,
286  std::map<int, std::vector<int>> const& ShowersMothers,
288  int planeid) const
289 {
290 
293 
294  //Find the energy for each track ID.
295  std::map<int, double> trackIDToEDepMap;
296  std::map<int, double> trackIDTo3EDepMap;
297  for (std::vector<art::Ptr<recob::Hit>>::const_iterator hitIt = hits.begin(); hitIt != hits.end();
298  ++hitIt) {
299  art::Ptr<recob::Hit> hit = *hitIt;
300 
301  //Get the plane ID
302  geo::WireID wireid = (*hitIt)->WireID();
303  int PlaneID = wireid.Plane;
304  std::vector<sim::TrackIDE> trackIDs = bt_serv->HitToTrackIDEs(clockData, hit);
305  for (unsigned int idIt = 0; idIt < trackIDs.size(); ++idIt) {
306  trackIDTo3EDepMap[std::abs(trackIDs[idIt].trackID)] += trackIDs[idIt].energy;
307  if (PlaneID == planeid) {
308  trackIDToEDepMap[std::abs(trackIDs[idIt].trackID)] += trackIDs[idIt].energy;
309  }
310  }
311  }
312 
313  //Find the energy for each showermother.
314  std::map<int, double> MotherIDtoEMap;
315  std::map<int, double> MotherIDto3EMap;
316  for (std::map<int, std::vector<int>>::const_iterator showermother = ShowersMothers.begin();
317  showermother != ShowersMothers.end();
318  ++showermother) {
319  for (std::vector<int>::const_iterator daughter = (showermother->second).begin();
320  daughter != (showermother->second).end();
321  ++daughter) {
322  MotherIDtoEMap[showermother->first] += trackIDToEDepMap[*daughter];
323  MotherIDto3EMap[showermother->first] += trackIDTo3EDepMap[*daughter];
324  }
325  }
326 
327  //Loop over the mothers to find the most like candidate by identifying which true shower deposited the most energy in the hits.
328  double maxenergy = -1;
329  int objectTrack = -99999;
330  for (std::map<int, double>::iterator mapIt = MotherIDto3EMap.begin();
331  mapIt != MotherIDto3EMap.end();
332  mapIt++) {
333  double energy_three = mapIt->second;
334  double trackid = mapIt->first;
335  if (energy_three > maxenergy) {
336  maxenergy = energy_three;
337  objectTrack = trackid;
338  }
339  }
340 
341  //If the none of the shower mother deposited no energy then we cannot match this.
342  if (maxenergy == 0) { return std::make_pair(-99999, -99999); }
343 
344  return std::make_pair(objectTrack, MotherIDtoEMap[objectTrack]);
345 }
Float_t x
Definition: compare.C:6
intermediate_table::iterator iterator
SubRunNumber_t subRun() const
Definition: Event.cc:35
int PdgCode() const
Definition: MCParticle.h:213
shower::LArPandoraShowerAlg fLArPandoraShowerAlg
size_t Self() const
Returns the index of this particle.
Definition: PFParticle.h:88
ROOT::Math::DisplacementVector3D< ROOT::Math::Cartesian3D< double >, ROOT::Math::GlobalCoordinateSystemTag > Vector_t
Type for representation of momenta in 3D space.
Definition: geo_vectors.h:160
std::vector< sim::TrackIDE > HitToTrackIDEs(detinfo::DetectorClocksData const &clockData, recob::Hit const &hit) const
int Mother() const
Definition: MCParticle.h:214
Float_t y
Definition: compare.C:6
std::map< int, const simb::MCParticle * > GetTrueParticleMap() const
list_type::const_iterator const_iterator
Definition: ParticleList.h:132
Double_t z
Definition: plot.C:276
constexpr auto abs(T v)
Returns the absolute value of the argument.
STL namespace.
intermediate_table::const_iterator const_iterator
LArPandoraShowerCheatingAlg(const fhicl::ParameterSet &pset)
MaybeLogger_< ELseverityLevel::ELsev_error, false > LogError
int TrackId() const
Definition: MCParticle.h:211
IDparameter< geo::PlaneID > PlaneID
Member type of validated geo::PlaneID parameter.
double x_min
Definition: berger.C:15
decltype(auto) constexpr end(T &&obj)
ADL-aware version of std::end.
Definition: StdUtils.h:77
auto vector(Vector const &v)
Returns a manipulator which will print the specified array.
Definition: DumpUtils.h:289
void hits()
Definition: readHits.C:15
art::ServiceHandle< cheat::ParticleInventoryService > particleInventory
parameter set interface
std::map< int, std::vector< int > > GetTrueChain(std::map< int, const simb::MCParticle * > &trueParticles) const
iterator begin()
Definition: ParticleList.h:305
bool CheckElement(const std::string &Name) const
void CheatDebugEVD(detinfo::DetectorClocksData const &clockData, const simb::MCParticle *trueParticle, art::Event const &Event, reco::shower::ShowerElementHolder &ShowerEleHolder, const art::Ptr< recob::PFParticle > &pfparticle) const
EventNumber_t event() const
Definition: Event.cc:41
PlaneID_t Plane
Index of the plane within its TPC.
Definition: geo_types.h:481
int GetElement(const std::string &Name, T &Element) const
int TrueParticleID(detinfo::DetectorClocksData const &clockData, const art::Ptr< recob::Hit > &hit) const
double SpacePointProjection(art::Ptr< recob::SpacePoint > const &sp, geo::Point_t const &vertex, geo::Vector_t const &direction) const
Detector simulation of raw signals on wires.
const sim::ParticleList & ParticleList() const
double x_max
Definition: berger.C:16
std::pair< int, double > TrueParticleIDFromTrueChain(detinfo::DetectorClocksData const &clockData, std::map< int, std::vector< int >> const &ShowersMothers, std::vector< art::Ptr< recob::Hit >> const &hits, int planeid) const
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:180
art::ServiceHandle< art::TFileService > tfs
ValidHandle< PROD > getValidHandle(InputTag const &tag) const
decltype(auto) get(T &&obj)
ADL-aware version of std::to_string.
Definition: StdUtils.h:120
Contains all timing reference information for the detector.
Float_t proj
Definition: plot.C:35
Definition: MVAAlg.h:12
constexpr WireID()=default
Default constructor: an invalid TPC ID.
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
Definition: Ptr.h:306
RunNumber_t run() const
Definition: Event.cc:29
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
Event finding and building.