#include <LArPandoraShowerCheatingAlg.h>

Public Member Functions
	LArPandoraShowerCheatingAlg (const fhicl::ParameterSet &pset)

std::map< int, const simb::MCParticle * >	GetTrueParticleMap () const

std::map< int, std::vector< int > >	GetTrueChain (std::map< int, const simb::MCParticle * > &trueParticles) 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

int	TrueParticleID (detinfo::DetectorClocksData const &clockData, const art::Ptr< recob::Hit > &hit) const

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

Private Attributes
shower::LArPandoraShowerAlg	fLArPandoraShowerAlg

art::InputTag	fHitModuleLabel

art::InputTag	fPFParticleLabel

art::ServiceHandle < cheat::ParticleInventoryService >	particleInventory

art::ServiceHandle < art::TFileService >	tfs

std::string	fShowerStartPositionInputLabel

std::string	fShowerDirectionInputLabel

std::string	fInitialTrackSpacePointsInputLabel

Detailed Description

Definition at line 33 of file LArPandoraShowerCheatingAlg.h.

Constructor & Destructor Documentation

shower::LArPandoraShowerCheatingAlg::LArPandoraShowerCheatingAlg ( const fhicl::ParameterSet & pset )

Definition at line 17 of file LArPandoraShowerCheatingAlg.cxx.

   : fLArPandoraShowerAlg(pset.get<fhicl::ParameterSet>("LArPandoraShowerAlg"))
   , fHitModuleLabel(pset.get<art::InputTag>("HitModuleLabel"))
   , fPFParticleLabel(pset.get<art::InputTag>("PFParticleLabel"))
   , fShowerStartPositionInputLabel(pset.get<std::string>("ShowerStartPositionInputLabel"))
   , fShowerDirectionInputLabel(pset.get<std::string>("ShowerDirectionInputLabel"))
   , fInitialTrackSpacePointsInputLabel(pset.get<std::string>("InitialTrackSpacePointsInputLabel"))
 {}

Member Function Documentation

void shower::LArPandoraShowerCheatingAlg::CheatDebugEVD	(	detinfo::DetectorClocksData const &	clockData,
		const simb::MCParticle *	trueParticle,
		art::Event const &	Event,
		reco::shower::ShowerElementHolder &	ShowerEleHolder,
		const art::Ptr< recob::PFParticle > &	pfparticle
	)		const

Definition at line 76 of file LArPandoraShowerCheatingAlg.cxx.

 {
 
   std::cout << "Making Debug Event Display" << std::endl;
 
   //Function for drawing reco showers to check direction and initial track selection
 
   // Get run info to make unique canvas names
   int run = Event.run();
   int subRun = Event.subRun();
   int event = Event.event();
   int PFPID = pfparticle->Self();
 
   // Create the canvas
   TString canvasName = Form("canvas_%i_%i_%i_%i", run, subRun, event, PFPID);
   TCanvas* canvas = tfs->make<TCanvas>(canvasName, canvasName);
 
   // Initialise variables
   double x = 0;
   double y = 0;
   double z = 0;
 
   std::vector<art::Ptr<recob::SpacePoint>> showerSpacePoints;
   std::vector<art::Ptr<recob::SpacePoint>> otherSpacePoints;
 
   auto const hitHandle = Event.getValidHandle<std::vector<recob::Hit>>(fHitModuleLabel);
   std::vector<art::Ptr<recob::Hit>> hits;
   art::fill_ptr_vector(hits, hitHandle);
 
   // Get the hits associated with the space points
   const art::FindManyP<recob::SpacePoint> fmsph(hitHandle, Event, fPFParticleLabel);
   if (!fmsph.isValid()) {
     throw cet::exception("LArPandoraShowerCheatingAlg")
       << "Spacepoint and hit association not valid. Stopping.";
   }
 
   std::map<art::Ptr<recob::SpacePoint>, art::Ptr<recob::Hit>> spacePointHitMap;
   //Get the hits from the true particle
   for (auto hit : hits) {
     int trueParticleID = std::abs(TrueParticleID(clockData, hit));
     std::vector<art::Ptr<recob::SpacePoint>> sps = fmsph.at(hit.key());
     if (sps.size() == 1) {
       art::Ptr<recob::SpacePoint> sp = sps.front();
       if (trueParticleID == trueParticle->TrackId()) { showerSpacePoints.push_back(sp); }
       else {
         otherSpacePoints.push_back(sp);
       }
     }
   }
 
   if (!ShowerEleHolder.CheckElement(fShowerStartPositionInputLabel)) {
     mf::LogError("LArPandoraShowerCheatingAlg")
       << "Start position not set, returning " << std::endl;
     return;
   }
   if (!ShowerEleHolder.CheckElement(fShowerDirectionInputLabel)) {
     mf::LogError("LArPandoraShowerCheatingAlg") << "Direction not set, returning " << std::endl;
     return;
   }
   if (!ShowerEleHolder.CheckElement(fInitialTrackSpacePointsInputLabel)) {
     mf::LogError("LArPandoraShowerCheatingAlg")
       << "TrackSpacePoints not set, returning " << std::endl;
     return;
   }
 
   // Get info from shower property holder
   TVector3 showerStartPosition = {-999, -999, -999};
   TVector3 showerDirection = {-999, -999, -999};
   std::vector<art::Ptr<recob::SpacePoint>> trackSpacePoints;
 
   ShowerEleHolder.GetElement(fShowerStartPositionInputLabel, showerStartPosition);
   ShowerEleHolder.GetElement(fShowerDirectionInputLabel, showerDirection);
   ShowerEleHolder.GetElement(fInitialTrackSpacePointsInputLabel, trackSpacePoints);
 
   // Create 3D point at vertex, chosed to be origin for ease of use of display
   double startXYZ[3] = {0, 0, 0};
   auto startPoly = std::make_unique<TPolyMarker3D>(1, startXYZ);
 
   // Get the min and max projections along the direction to know how long to draw
   // the direction line
   double minProj = std::numeric_limits<double>::max();
   double maxProj = -std::numeric_limits<double>::max();
 
   double x_min = std::numeric_limits<double>::max(), x_max = -std::numeric_limits<double>::max();
   double y_min = std::numeric_limits<double>::max(), y_max = -std::numeric_limits<double>::max();
   double z_min = std::numeric_limits<double>::max(), z_max = -std::numeric_limits<double>::max();
 
   //initialise counter point
   int point = 0;
 
   // Make 3D points for each spacepoint in the shower
   auto showerPoly = std::make_unique<TPolyMarker3D>(showerSpacePoints.size());
   for (auto spacePoint : showerSpacePoints) {
     TVector3 pos = fLArPandoraShowerAlg.SpacePointPosition(spacePoint) - showerStartPosition;
 
     x = pos.X();
     y = pos.Y();
     z = pos.Z();
     x_min = std::min(x, x_min);
     x_max = std::max(x, x_max);
     y_min = std::min(y, y_min);
     y_max = std::max(y, y_max);
     z_min = std::min(z, z_min);
     z_max = std::max(z, z_max);
 
     showerPoly->SetPoint(point, x, y, z);
     ++point;
 
     // Calculate the projection of (point-startpoint) along the direction
     double proj =
       fLArPandoraShowerAlg.SpacePointProjection(spacePoint, showerStartPosition, showerDirection);
 
     maxProj = std::max(proj, maxProj);
     minProj = std::min(proj, minProj);
   } // loop over spacepoints
 
   // Create TPolyLine3D arrays
   double xDirPoints[2] = {minProj * showerDirection.X(), maxProj * showerDirection.X()};
   double yDirPoints[2] = {minProj * showerDirection.Y(), maxProj * showerDirection.Y()};
   double zDirPoints[2] = {minProj * showerDirection.Z(), maxProj * showerDirection.Z()};
 
   auto dirPoly = std::make_unique<TPolyLine3D>(2, xDirPoints, yDirPoints, zDirPoints);
 
   point = 0; // re-initialise counter
   auto trackPoly = std::make_unique<TPolyMarker3D>(trackSpacePoints.size());
   for (auto spacePoint : trackSpacePoints) {
     TVector3 pos = fLArPandoraShowerAlg.SpacePointPosition(spacePoint) - showerStartPosition;
     x = pos.X();
     y = pos.Y();
     z = pos.Z();
     trackPoly->SetPoint(point, x, y, z);
     ++point;
   } // loop over track spacepoints
 
   //  we want to draw all of the PFParticles in the event
   //Get the PFParticles
 
   auto otherPoly = std::make_unique<TPolyMarker3D>(otherSpacePoints.size());
 
   // initialise counters
   point = 0; // re-initialise counter
 
   for (auto const& sp : otherSpacePoints) {
     TVector3 pos = fLArPandoraShowerAlg.SpacePointPosition(sp) - showerStartPosition;
     x = pos.X();
     y = pos.Y();
     z = pos.Z();
     x_min = std::min(x, x_min);
     x_max = std::max(x, x_max);
     y_min = std::min(y, y_min);
     y_max = std::max(y, y_max);
     z_min = std::min(z, z_min);
     z_max = std::max(z, z_max);
     otherPoly->SetPoint(point, x, y, z);
     ++point;
   }
 
   gStyle->SetOptStat(0);
   TH3F axes("axes", "", 1, x_min, x_max, 1, y_min, y_max, 1, z_min, z_max);
   axes.SetDirectory(0);
   axes.GetXaxis()->SetTitle("X");
   axes.GetYaxis()->SetTitle("Y");
   axes.GetZaxis()->SetTitle("Z");
   axes.Draw();
 
   otherPoly->SetMarkerStyle(20);
   otherPoly->SetMarkerColor(4);
   otherPoly->Draw();
 
   // Draw all of the things
   showerPoly->SetMarkerStyle(20);
   showerPoly->Draw();
   trackPoly->SetMarkerStyle(20);
   trackPoly->SetMarkerColor(2);
   trackPoly->Draw();
   startPoly->SetMarkerStyle(21);
   startPoly->SetMarkerSize(2);
   startPoly->SetMarkerColor(3);
   startPoly->Draw();
   dirPoly->SetLineWidth(1);
   dirPoly->SetLineColor(6);
   dirPoly->Draw();
 
   // Save the canvas. Don't usually need this when using TFileService but this in the alg
   // not a module and didn't work without this so im going with it.
   canvas->Write();
 }

std::map< int, std::vector< int > > shower::LArPandoraShowerCheatingAlg::GetTrueChain ( std::map< int, const simb::MCParticle * > & trueParticles ) const

Definition at line 45 of file LArPandoraShowerCheatingAlg.cxx.

 {
 
   // Roll up showers if not already done:
   std::map<int, std::vector<int>> showerMothers;
 
   // Loop over daughters and find th`e mothers
   for (const auto& particleIt : trueParticles) {
     const simb::MCParticle* particle = particleIt.second;
     const simb::MCParticle* mother = particle;
 
     if (std::abs(particle->PdgCode()) != 11 && std::abs(particle->PdgCode()) != 22) { continue; }
 
     // While the grand mother exists and is an electron or photon
     // Note the true mother will skip this loop and fill itself into the map
     while (mother->Mother() != 0 && trueParticles.find(mother->Mother()) != trueParticles.end()) {
 
       int motherId = mother->Mother();
       if (std::abs(trueParticles[motherId]->PdgCode()) != 11 &&
           std::abs(trueParticles[motherId]->PdgCode()) != 22) {
         break;
       }
       mother = trueParticles[motherId];
     }
     showerMothers[mother->TrackId()].push_back(particle->TrackId());
   }
   return showerMothers;
 }

std::map< int, const simb::MCParticle * > shower::LArPandoraShowerCheatingAlg::GetTrueParticleMap ( ) const

Definition at line 27 of file LArPandoraShowerCheatingAlg.cxx.

 {
 
   const sim::ParticleList& particles = particleInventory->ParticleList();
 
   std::map<int, const simb::MCParticle*> trueParticles;
   // Make a map of track id to particle
   for (sim::ParticleList::const_iterator particleIt = particles.begin();
        particleIt != particles.end();
        ++particleIt) {
     const simb::MCParticle* particle = particleIt->second;
     trueParticles[particle->TrackId()] = particle;
     //std::cout<<"Particle ID: "<<particle->TrackId()<<" and PDG: "<<particle->PdgCode()<<std::endl;
   }
   return trueParticles;
 }

int shower::LArPandoraShowerCheatingAlg::TrueParticleID	(	detinfo::DetectorClocksData const &	clockData,
		const art::Ptr< recob::Hit > &	hit
	)		const

Definition at line 270 of file LArPandoraShowerCheatingAlg.cxx.

 {
 
   double particleEnergy = 0;
   int likelyTrackID = 0;
   art::ServiceHandle<cheat::BackTrackerService> bt_serv;
   std::vector<sim::TrackIDE> trackIDs = bt_serv->HitToTrackIDEs(clockData, hit);
   for (unsigned int idIt = 0; idIt < trackIDs.size(); ++idIt) {
     if (trackIDs.at(idIt).energy > particleEnergy) {
       particleEnergy = trackIDs.at(idIt).energy;
       likelyTrackID = trackIDs.at(idIt).trackID;
     }
   }
   return likelyTrackID;
 }

std::pair< int, double > shower::LArPandoraShowerCheatingAlg::TrueParticleIDFromTrueChain	(	detinfo::DetectorClocksData const &	clockData,
		std::map< int, std::vector< int >> const &	ShowersMothers,
		std::vector< art::Ptr< recob::Hit >> const &	hits,
		int	planeid
	)		const

Definition at line 288 of file LArPandoraShowerCheatingAlg.cxx.

 {
 
   art::ServiceHandle<cheat::BackTrackerService> bt_serv;
   art::ServiceHandle<cheat::ParticleInventoryService> particleInventory;
 
   //Find the energy for each track ID.
   std::map<int, double> trackIDToEDepMap;
   std::map<int, double> trackIDTo3EDepMap;
   for (std::vector<art::Ptr<recob::Hit>>::const_iterator hitIt = hits.begin(); hitIt != hits.end();
        ++hitIt) {
     art::Ptr<recob::Hit> hit = *hitIt;
 
     //Get the plane ID
     geo::WireID wireid = (*hitIt)->WireID();
     int PlaneID = wireid.Plane;
     std::vector<sim::TrackIDE> trackIDs = bt_serv->HitToTrackIDEs(clockData, hit);
     for (unsigned int idIt = 0; idIt < trackIDs.size(); ++idIt) {
       trackIDTo3EDepMap[std::abs(trackIDs[idIt].trackID)] += trackIDs[idIt].energy;
       if (PlaneID == planeid) {
         trackIDToEDepMap[std::abs(trackIDs[idIt].trackID)] += trackIDs[idIt].energy;
       }
     }
   }
 
   //Find the energy for each showermother.
   std::map<int, double> MotherIDtoEMap;
   std::map<int, double> MotherIDto3EMap;
   for (std::map<int, std::vector<int>>::const_iterator showermother = ShowersMothers.begin();
        showermother != ShowersMothers.end();
        ++showermother) {
     for (std::vector<int>::const_iterator daughter = (showermother->second).begin();
          daughter != (showermother->second).end();
          ++daughter) {
       MotherIDtoEMap[showermother->first] += trackIDToEDepMap[*daughter];
       MotherIDto3EMap[showermother->first] += trackIDTo3EDepMap[*daughter];
     }
   }
 
   //Loop over the mothers to find the most like candidate by identifying which true shower deposited the most energy in the hits.
   double maxenergy = -1;
   int objectTrack = -99999;
   for (std::map<int, double>::iterator mapIt = MotherIDto3EMap.begin();
        mapIt != MotherIDto3EMap.end();
        mapIt++) {
     double energy_three = mapIt->second;
     double trackid = mapIt->first;
     if (energy_three > maxenergy) {
       maxenergy = energy_three;
       objectTrack = trackid;
     }
   }
 
   //If the none of the shower mother deposited no energy then we cannot match this.
   if (maxenergy == 0) { return std::make_pair(-99999, -99999); }
 
   return std::make_pair(objectTrack, MotherIDtoEMap[objectTrack]);
 }