10 #include "TPolyLine3D.h" 11 #include "TPolyMarker3D.h" 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"))
32 std::map<int, const simb::MCParticle*> trueParticles;
35 particleIt != particles.
end();
38 trueParticles[particle->
TrackId()] = particle;
45 std::map<int, const simb::MCParticle*>& trueParticles)
const 49 std::map<int, std::vector<int>> showerMothers;
52 for (
const auto& particleIt : trueParticles) {
60 while (mother->
Mother() != 0 && trueParticles.find(mother->
Mother()) != trueParticles.end()) {
62 int motherId = mother->
Mother();
63 if (
std::abs(trueParticles[motherId]->PdgCode()) != 11 &&
64 std::abs(trueParticles[motherId]->PdgCode()) != 22) {
67 mother = trueParticles[motherId];
82 std::cout <<
"Making Debug Event Display" << std::endl;
87 int run = Event.
run();
88 int subRun = Event.
subRun();
89 int event = Event.
event();
90 int PFPID = pfparticle->
Self();
93 TString canvasName = Form(
"canvas_%i_%i_%i_%i", run, subRun,
event, PFPID);
94 TCanvas* canvas =
tfs->make<TCanvas>(canvasName, canvasName);
101 std::vector<art::Ptr<recob::SpacePoint>> showerSpacePoints;
102 std::vector<art::Ptr<recob::SpacePoint>> otherSpacePoints;
105 std::vector<art::Ptr<recob::Hit>>
hits;
110 if (!fmsph.isValid()) {
112 <<
"Spacepoint and hit association not valid. Stopping.";
117 for (
auto hit : hits) {
119 std::vector<art::Ptr<recob::SpacePoint>> sps = fmsph.at(
hit.key());
120 if (sps.size() == 1) {
122 if (trueParticleID == trueParticle->
TrackId()) { showerSpacePoints.push_back(sp); }
124 otherSpacePoints.push_back(sp);
131 <<
"Start position not set, returning " << std::endl;
135 mf::LogError(
"LArPandoraShowerCheatingAlg") <<
"Direction not set, returning " << std::endl;
140 <<
"TrackSpacePoints not set, returning " << std::endl;
147 std::vector<art::Ptr<recob::SpacePoint>> trackSpacePoints;
154 double startXYZ[3] = {0, 0, 0};
155 auto startPoly = std::make_unique<TPolyMarker3D>(1, startXYZ);
159 double minProj = std::numeric_limits<double>::max();
160 double maxProj = -std::numeric_limits<double>::max();
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();
170 auto showerPoly = std::make_unique<TPolyMarker3D>(showerSpacePoints.size());
171 for (
auto const& spacePoint : showerSpacePoints) {
172 auto const pos = spacePoint->position() - showerStartPosition;
176 x_min = std::min(x, x_min);
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);
183 showerPoly->SetPoint(point, x, y, z);
190 maxProj = std::max(proj, maxProj);
191 minProj = std::min(proj, minProj);
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()};
199 auto dirPoly = std::make_unique<TPolyLine3D>(2, xDirPoints, yDirPoints, zDirPoints);
202 auto trackPoly = std::make_unique<TPolyMarker3D>(trackSpacePoints.size());
203 for (
auto const& spacePoint : trackSpacePoints) {
204 auto const pos = spacePoint->position() - showerStartPosition;
208 trackPoly->SetPoint(point, x, y, z);
215 auto otherPoly = std::make_unique<TPolyMarker3D>(otherSpacePoints.size());
220 for (
auto const& sp : otherSpacePoints) {
221 auto const pos = sp->position() - showerStartPosition;
225 x_min = std::min(x, x_min);
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);
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");
243 otherPoly->SetMarkerStyle(20);
244 otherPoly->SetMarkerColor(4);
248 showerPoly->SetMarkerStyle(20);
250 trackPoly->SetMarkerStyle(20);
251 trackPoly->SetMarkerColor(2);
253 startPoly->SetMarkerStyle(21);
254 startPoly->SetMarkerSize(2);
255 startPoly->SetMarkerColor(3);
257 dirPoly->SetLineWidth(1);
258 dirPoly->SetLineColor(6);
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;
281 return likelyTrackID;
286 std::map<
int, std::vector<int>>
const& ShowersMothers,
295 std::map<int, double> trackIDToEDepMap;
296 std::map<int, double> trackIDTo3EDepMap;
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;
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();
320 daughter != (showermother->second).
end();
322 MotherIDtoEMap[showermother->first] += trackIDToEDepMap[*daughter];
323 MotherIDto3EMap[showermother->first] += trackIDTo3EDepMap[*daughter];
328 double maxenergy = -1;
329 int objectTrack = -99999;
331 mapIt != MotherIDto3EMap.end();
333 double energy_three = mapIt->second;
334 double trackid = mapIt->first;
335 if (energy_three > maxenergy) {
336 maxenergy = energy_three;
337 objectTrack = trackid;
342 if (maxenergy == 0) {
return std::make_pair(-99999, -99999); }
344 return std::make_pair(objectTrack, MotherIDtoEMap[objectTrack]);
SubRunNumber_t subRun() const
shower::LArPandoraShowerAlg fLArPandoraShowerAlg
size_t Self() const
Returns the index of this particle.
ROOT::Math::DisplacementVector3D< ROOT::Math::Cartesian3D< double >, ROOT::Math::GlobalCoordinateSystemTag > Vector_t
Type for representation of momenta in 3D space.
std::vector< sim::TrackIDE > HitToTrackIDEs(detinfo::DetectorClocksData const &clockData, recob::Hit const &hit) const
std::map< int, const simb::MCParticle * > GetTrueParticleMap() const
list_type::const_iterator const_iterator
constexpr auto abs(T v)
Returns the absolute value of the argument.
LArPandoraShowerCheatingAlg(const fhicl::ParameterSet &pset)
MaybeLogger_< ELseverityLevel::ELsev_error, false > LogError
IDparameter< geo::PlaneID > PlaneID
Member type of validated geo::PlaneID parameter.
decltype(auto) constexpr end(T &&obj)
ADL-aware version of std::end.
auto vector(Vector const &v)
Returns a manipulator which will print the specified array.
art::ServiceHandle< cheat::ParticleInventoryService > particleInventory
std::map< int, std::vector< int > > GetTrueChain(std::map< int, const simb::MCParticle * > &trueParticles) const
std::string fShowerDirectionInputLabel
bool CheckElement(const std::string &Name) const
art::InputTag fHitModuleLabel
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
PlaneID_t Plane
Index of the plane within its TPC.
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
std::string fShowerStartPositionInputLabel
Detector simulation of raw signals on wires.
const sim::ParticleList & ParticleList() const
std::string fInitialTrackSpacePointsInputLabel
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.
art::ServiceHandle< art::TFileService > tfs
ValidHandle< PROD > getValidHandle(InputTag const &tag) const
decltype(auto) get(T &&obj)
ADL-aware version of std::to_string.
Contains all timing reference information for the detector.
constexpr WireID()=default
Default constructor: an invalid TPC ID.
art::InputTag fPFParticleLabel
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
cet::coded_exception< error, detail::translate > exception
Event finding and building.