pizero::ShowerProcess Class Reference

#include <ShowerProcess.h>

Public Member Functions
	ShowerProcess (const simb::MCParticle &mcpart, const art::Event &evt, detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, std::string showerLabel="pandoraShower", std::string trackLabel="pandoraTrack", std::string simulationLabel="largeant")
	ShowerProcess (const recob::Shower &shower, const art::Event &evt, detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, std::string showerLabel="pandoraShower", std::string trackLabel="pandoraTrack", std::string simulationLabel="largeant")
	ShowerProcess (const simb::MCParticle &mcpart, const recob::Shower &shower, const art::Event &evt, detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, std::string showerLabel="pandoraShower", std::string trackLabel="pandoraTrack", std::string simulationLabel="largeant")
	~ShowerProcess ()
const simb::MCParticle *	mcparticle () const
const recob::Shower *	shower () const
const recob::Track *	track () const
const Cone *	cone () const
std::vector< const simb::MCParticle * >	mcparticles () const
std::vector< const recob::Shower * >	showers () const
std::vector< const recob::Track * >	tracks () const
double	energy (const std::string &caloLabel, const double calib, const double norm, const std::string &scefile, const std::string &yzfile, const std::string &xfile) const

Private Member Functions
void	find_mcparticles ()
void	find_showers ()
void	find_tracks ()
void	fill_cone ()

Private Attributes
std::vector< const simb::MCParticle * >	m_mcparts
std::vector< const recob::Shower * >	m_showers
std::vector< const recob::Track * >	m_tracks
Cone *	m_cone = nullptr
const art::Event *	m_evt
detinfo::DetectorClocksData const *	m_clockdata
detinfo::DetectorPropertiesData const *	m_detprop
std::string	m_showerLabel
std::string	m_trackLabel
std::string	m_simulationLabel
protoana::ProtoDUNETruthUtils	truthUtils
protoana::ProtoDUNEShowerUtils	shUtils

Detailed Description

Definition at line 25 of file ShowerProcess.h.

Constructor & Destructor Documentation

pizero::ShowerProcess::ShowerProcess ( const simb::MCParticle &  mcpart, const art::Event &  evt, detinfo::DetectorClocksData const &  clockData, detinfo::DetectorPropertiesData const &  detProp, std::string  showerLabel = "pandoraShower", std::string  trackLabel = "pandoraTrack", std::string  simulationLabel = "largeant"  )

inline

Definition at line 55 of file ShowerProcess.h.

                                                       :
                m_evt(&evt),
                m_clockdata(&clockData),
                m_detprop(&detProp),
                m_showerLabel(showerLabel),
                m_trackLabel(trackLabel),
                m_simulationLabel(simulationLabel) {
    m_mcparts.push_back(&mcpart);
    find_showers();
    find_tracks();
    fill_cone();
  }

pizero::ShowerProcess::ShowerProcess ( const recob::Shower &  shower, const art::Event &  evt, detinfo::DetectorClocksData const &  clockData, detinfo::DetectorPropertiesData const &  detProp, std::string  showerLabel = "pandoraShower", std::string  trackLabel = "pandoraTrack", std::string  simulationLabel = "largeant"  )

inline

Definition at line 73 of file ShowerProcess.h.

                                                       :
                m_evt(&evt),
                m_clockdata(&clockData),
                m_detprop(&detProp),
                m_showerLabel(showerLabel),
                m_trackLabel(trackLabel),
                m_simulationLabel(simulationLabel) {
    m_showers.push_back(&shower);
    find_mcparticles();
    find_tracks();
    fill_cone();
  }

pizero::ShowerProcess::ShowerProcess ( const simb::MCParticle &  mcpart, const recob::Shower &  shower, const art::Event &  evt, detinfo::DetectorClocksData const &  clockData, detinfo::DetectorPropertiesData const &  detProp, std::string  showerLabel = "pandoraShower", std::string  trackLabel = "pandoraTrack", std::string  simulationLabel = "largeant"  )

inline

Definition at line 90 of file ShowerProcess.h.

                                                       :
                m_evt(&evt),
                m_clockdata(&clockData),
                m_detprop(&detProp),
                m_showerLabel(showerLabel),
                m_trackLabel(trackLabel),
                m_simulationLabel(simulationLabel) {
    m_mcparts.push_back(&mcpart);
    m_showers.push_back(&shower);
    find_tracks();
    fill_cone();
  }

pizero::ShowerProcess::~ShowerProcess ( )

inline

Definition at line 108 of file ShowerProcess.h.

{ if(m_cone != nullptr) delete m_cone; } // TODO: make unique_ptr.

Member Function Documentation

const Cone* pizero::ShowerProcess::cone ( ) const

inline

Definition at line 120 of file ShowerProcess.h.

                           {
    return m_cone;
  }

double pizero::ShowerProcess::energy	(	const std::string &	caloLabel,
		const double	calib,
		const double	norm,
		const std::string &	scefile,
		const std::string &	yzfile,
		const std::string &	xfile
	)		const

Definition at line 260 of file ShowerProcess.h.

                                {
  double totE = 0;

  // Check whether there is a shower associated with this object.
  if(shower() == nullptr) return totE;

  // Constants from DUNE docDB 15974 by A Paudel.
  const double Wion = 23.6e-6; // ArgoNeuT-determined parameter at E0 kV/cm
  const double recob = 0.6417; // Recombination factor from Aaron's talk.

  // Get dQ/dx YZ and X correction factors.
  TFile* yzf = new TFile(yzfile.c_str());
  TH2F* yzcorrposhist =(TH2F*)yzf->Get("correction_dqdx_ZvsY_positiveX_hist_2");
  TH2F* yzcorrneghist =(TH2F*)yzf->Get("correction_dqdx_ZvsY_negativeX_hist_2");
  TFile* xf = new TFile(xfile.c_str());
  TH1F* xcorrhist = (TH1F*)xf->Get("dqdx_X_correction_hist_2");

  // Get the calorimetry objects from the event using the shower utilities.
  std::vector<anab::Calorimetry> calovec =
    shUtils.GetRecoShowerCalorimetry(*shower(), *m_evt, m_showerLabel, caloLabel);

  // Loop over all plane calorimetry objects, but only use the collection plane.
  for(const anab::Calorimetry& calo : calovec) {
    if(calo.PlaneID().Plane != 2) continue;

    // Loop over all hits in the collection plane calorimetry object.
    for(unsigned i = 0; i < calo.dQdx().size(); ++i) {

      double dQdx = calo.dQdx()[i];
      const double pitch = calo.TrkPitchVec()[i];
      const geo::Point_t& pos = calo.XYZ()[i];
      const double x = pos.X();
      const double y = pos.Y();
      const double z = pos.Z();
      // Can't apply location-dependent corrections if there's no location available.
      if(x==0 && y==0 && z==0) {
        totE += dQdx*pitch*norm/calib*Wion/recob * 1e-3;
        continue;
      }
      // Spatial corrections.
      double yzcorr;
      if(x >= 0) {
        yzcorr = yzcorrposhist->GetBinContent(yzcorrposhist->FindBin(z,y));
      } else {
        yzcorr = yzcorrneghist->GetBinContent(yzcorrneghist->FindBin(z,y));
      }
      // Correct dQ/dx.
      // std::cout << "Bare energy (GeV): " << dQdx*pitch*norm/calib*Wion/recob * 1e-3 << '\n';
      const double xcorr = xcorrhist->GetBinContent(xcorrhist->FindBin(x));
      // std::cout << "Hit energy (GeV): " << xcorr*yzcorr*dQdx*pitch*norm/calib*Wion/recob * 1e-3 << "\n\n";
      totE += xcorr*yzcorr*dQdx*pitch*norm/calib*Wion/recob * 1e-3; // Convert from MeV to GeV.
    } // for hits in calo
  } // for calo objects

  // std::cout << "Total energy: " << totE << '\n';
  return totE;
}

void pizero::ShowerProcess::fill_cone ( )

private

Definition at line 242 of file ShowerProcess.h.

                              {
  // TODO: Find the start and direction of first reconstructed object.
  TVector3 cstart, cdir;
  if(m_showers.size() != 0) {
    cstart = m_showers[0]->ShowerStart();
    cdir = m_showers[0]->Direction();
    // Make cone object
    m_cone = new Cone(*m_evt, *m_clockdata, *m_detprop, cstart, cdir, m_showers[0]->Length());
  }
  // else if(m_tracks.size() != 0) {
  //   const recob::tracking::Point_t tstart = m_tracks[0]->Start();
  //   const recob::tracking::Vector_t tdir = m_tracks[0]->StartDirection();
  //   cstart = TVector3(tstart.X(), tstart.Y(), tstart.Z());
  //   cdir = TVector3(tdir.X(), tdir.Y(), tdir.Z());
  // }

}

void pizero::ShowerProcess::find_mcparticles ( )

private

Definition at line 135 of file ShowerProcess.h.

                                     {
  if(m_evt->isRealData() || m_showers.size() == 0) return;
  m_mcparts.push_back(truthUtils.GetMCParticleFromReco(*m_clockdata, *m_showers[0], *m_evt, m_showerLabel));
}

void pizero::ShowerProcess::find_showers ( )

private

Definition at line 141 of file ShowerProcess.h.

                                 {
  if(m_mcparts.size() == 0) return;

  // Find all showers this MCParticle contributed to.
  std::vector<std::pair<const recob::Shower*, double>> showers =
    truthUtils.GetRecoShowerListFromMCParticle(*m_clockdata, *m_mcparts[0], *m_evt, m_showerLabel);
  // Only showers to which the MCParticle was the primary contributor are saved.
  for(const std::pair<const recob::Shower*, double>& sh : showers) {
    const simb::MCParticle* sh_part =
      truthUtils.GetMCParticleFromReco(*m_clockdata, *sh.first, *m_evt, m_showerLabel);
    if(std::abs(sh_part->TrackId()) == std::abs(m_mcparts[0]->TrackId())) {
      m_showers.push_back(sh.first);
    }
  }
}

void pizero::ShowerProcess::find_tracks ( )

private

Definition at line 158 of file ShowerProcess.h.

                                {
  // For now return if no MC present.
  if(m_mcparts.size() == 0) return;

  using weightedTrackPair = std::pair<const recob::Track*, double>;
  std::vector<weightedTrackPair> outVec;

  // Get all reconstructed tracks
  auto allRecoTracks = m_evt->getValidHandle<std::vector<recob::Track> >(m_trackLabel);

  // We need the association between the tracks and the hits
  const art::FindManyP<recob::Hit> findTrackHits(allRecoTracks, *m_evt, m_trackLabel);

  art::ServiceHandle<cheat::BackTrackerService> bt_serv;
  art::ServiceHandle<cheat::ParticleInventoryService> pi_serv;
  std::unordered_map<int, double> recoTrack;

  // Record the energy contribution to the MCParticle of every relevant reco track
  double part_E = 0;
  for (recob::Track const & track : *allRecoTracks) {
    // Record the contribution of all MCParticles to this track.
    std::unordered_map<int, double> mcpmap;
    // Loop through the track's hits.
    for (art::Ptr<recob::Hit> const & hptr : findTrackHits.at(track.ID())) {
      // Loop over hit IDEs to find our MCParticle's part in it
      for (const sim::IDE* ide : bt_serv->HitToSimIDEs_Ps(*m_clockdata, hptr)) {
        mcpmap[abs(ide->trackID)] += ide->energy;
        if (abs(ide->trackID) == mcparticle()->TrackId()) {
          recoTrack[track.ID()] += ide->energy;
        }
        part_E += ide->energy;
      } // sim::IDE*
    } // art::Ptr<recob::Hit>

    // Remove track if the shower's MCParticle was not the main contributor.
    if(recoTrack.size() > 0) {
      int maxMCPID = mcpmap[0];
      if(mcpmap.size() > 1) {
        maxMCPID = std::max_element(mcpmap.begin(), mcpmap.end(),
          [](const std::pair<int, double>& a, const std::pair<int, double>& b){
            return a.second < b.second;
          })->first;
      }
      if(maxMCPID != mcparticle()->TrackId()) {
        recoTrack.erase(track.ID());
      }
    }
  } // const recob::Track&
  // for(auto const& p : recoTrack) {
  //   std::cout << "Track with ID " << p.first << " has energy " << p.second << '\n';
  // }

  // Fill and sort the output vector
  for (std::pair<int, double> const& p : recoTrack)
  {
    auto const trackIt = std::find_if(allRecoTracks->begin(), allRecoTracks->end(),
                               [&](recob::Track tr){ return tr.ID() == p.first; });
    outVec.push_back(std::make_pair(&*trackIt, p.second));
  }
  std::sort(outVec.begin(), outVec.end(),
    [](weightedTrackPair a, weightedTrackPair b){ return a.second > b.second;});

  // Normalise the vector weights
  if (part_E < 1e-5) { part_E = 1; } // Protect against zero division
  for (weightedTrackPair& p : outVec)
  {
    p.second /= part_E;
    m_tracks.push_back(p.first);
  }

  // // Find all tracks this MCParticle contributed to.
  // std::vector<std::pair<const recob::Track*, double>> tracks =
  //   truthUtils.GetRecoTrackListFromMCParticle(*m_mcparts[0], *m_evt, m_trackLabel);
  // // Only tracks to which the MCParticle was the primary contributor are saved.
  // for(const std::pair<const recob::Track*, double>& tr : tracks) {
  //   const simb::MCParticle* tr_part =
  //     truthUtils.GetMCParticleFromReco(*tr.first, *m_evt, m_trackLabel);
  //   if(std::abs(tr_part->TrackId()) == std::abs(m_mcparts[0]->TrackId())) {
  //     m_tracks.push_back(tr.first);
  //   }
  // }
}

const simb::MCParticle* pizero::ShowerProcess::mcparticle ( ) const

inline

Definition at line 111 of file ShowerProcess.h.

                                           {
    return m_mcparts.size()!=0? m_mcparts[0]: nullptr;
  }

std::vector<const simb::MCParticle*> pizero::ShowerProcess::mcparticles ( ) const

inline

Definition at line 124 of file ShowerProcess.h.

{ return m_mcparts; }

const recob::Shower* pizero::ShowerProcess::shower ( ) const

inline

Definition at line 114 of file ShowerProcess.h.

                                    {
    return m_showers.size()!=0? m_showers[0]: nullptr;
  }

std::vector<const recob::Shower*> pizero::ShowerProcess::showers ( ) const

inline

Definition at line 125 of file ShowerProcess.h.

{ return m_showers; }

const recob::Track* pizero::ShowerProcess::track ( ) const

inline

Definition at line 117 of file ShowerProcess.h.

                                  {
    return m_tracks.size()!=0? m_tracks[0]: nullptr;
  }

std::vector<const recob::Track*> pizero::ShowerProcess::tracks ( ) const

inline

Definition at line 126 of file ShowerProcess.h.

{ return m_tracks; }

Member Data Documentation

detinfo::DetectorClocksData const* pizero::ShowerProcess::m_clockdata

private

Definition at line 37 of file ShowerProcess.h.

Cone* pizero::ShowerProcess::m_cone = nullptr

private

Definition at line 33 of file ShowerProcess.h.

detinfo::DetectorPropertiesData const* pizero::ShowerProcess::m_detprop

private

Definition at line 38 of file ShowerProcess.h.

const art::Event* pizero::ShowerProcess::m_evt

private

Definition at line 36 of file ShowerProcess.h.

std::vector<const simb::MCParticle*> pizero::ShowerProcess::m_mcparts

private

Definition at line 28 of file ShowerProcess.h.

std::string pizero::ShowerProcess::m_showerLabel

private

Definition at line 39 of file ShowerProcess.h.

std::vector<const recob::Shower*> pizero::ShowerProcess::m_showers

private

Definition at line 29 of file ShowerProcess.h.

std::string pizero::ShowerProcess::m_simulationLabel

private

Definition at line 41 of file ShowerProcess.h.

std::string pizero::ShowerProcess::m_trackLabel

private

Definition at line 40 of file ShowerProcess.h.

std::vector<const recob::Track*> pizero::ShowerProcess::m_tracks

private

Definition at line 30 of file ShowerProcess.h.

protoana::ProtoDUNEShowerUtils pizero::ShowerProcess::shUtils

private

Definition at line 51 of file ShowerProcess.h.

protoana::ProtoDUNETruthUtils pizero::ShowerProcess::truthUtils

private

Definition at line 50 of file ShowerProcess.h.

---

The documentation for this class was generated from the following file:

• protoduneana/protoduneana/singlephase/EMTaskForce/ShowerProcess.h