#include <PointIdAlg.h>

Inheritance diagram for nnet::TrainingDataAlg:

Classes
struct	Config

struct	WireDrift

Public Types
enum	EMask { kNone = 0, kPdgMask = 0x00000FFF, kTypeMask = 0x0000F000, kVtxMask = 0xFFFF0000 }

enum	ETrkType { kDelta = 0x1000, kMichel = 0x2000, kPriEl = 0x4000, kPriMu = 0x8000 }

enum	EVtxId { kNuNC = 0x0010000, kNuCC = 0x0020000, kNuPri = 0x0040000, kNuE = 0x0100000, kNuMu = 0x0200000, kNuTau = 0x0400000, kHadr = 0x1000000, kPi0 = 0x2000000, kDecay = 0x4000000, kConv = 0x8000000, kElectronEnd = 0x10000000, kElastic = 0x20000000, kInelastic = 0x40000000 }

Public Types inherited from img::DataProviderAlg
enum	EDownscaleMode { kMax = 1, kMaxMean = 2, kMean = 3 }

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

	TrainingDataAlg (const Config &config)

	~TrainingDataAlg () override

void	reconfigure (const Config &config)

bool	saveSimInfo () const

bool	setEventData (const art::Event &event, detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, unsigned int plane, unsigned int tpc, unsigned int cryo)

bool	setDataEventData (const art::Event &event, detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, unsigned int plane, unsigned int tpc, unsigned int cryo)

bool	findCrop (float max_e_cut, unsigned int &w0, unsigned int &w1, unsigned int &d0, unsigned int &d1) const

double	getEdepTot () const

std::vector< float > const &	wireEdep (size_t widx) const

std::vector< int > const &	wirePdg (size_t widx) const

Public Member Functions inherited from img::DataProviderAlg
	DataProviderAlg (const fhicl::ParameterSet &pset)

	DataProviderAlg (const Config &config)

virtual	~DataProviderAlg ()

bool	setWireDriftData (const detinfo::DetectorClocksData &clock_data, const detinfo::DetectorPropertiesData &det_prop, const std::vector< recob::Wire > &wires, unsigned int plane, unsigned int tpc, unsigned int cryo)

std::vector< float > const &	wireData (size_t widx) const

std::vector< std::vector< float > >	getPatch (size_t wire, float drift, size_t patchSizeW, size_t patchSizeD) const

float	getPixelOrZero (int wire, int drift) const

double	getAdcSum () const

size_t	getAdcArea () const

float	poolMax (int wire, int drift, size_t r=0) const
	Pool max value in a patch around the wire/drift pixel. More...

unsigned int	Cryo () const
	Pool sum of pixels in a patch around the wire/drift pixel. More...

unsigned int	TPC () const

unsigned int	Plane () const

unsigned int	NWires () const

unsigned int	NScaledDrifts () const

unsigned int	NCachedDrifts () const

unsigned int	DriftWindow () const

float	ZeroLevel () const
	Level of zero ADC after scaling. More...

double	LifetimeCorrection (detinfo::DetectorClocksData const &clock_data, detinfo::DetectorPropertiesData const &det_prop, double tick) const

Protected Member Functions
img::DataProviderAlgView	resizeView (detinfo::DetectorClocksData const &clock_data, detinfo::DetectorPropertiesData const &det_prop, size_t wires, size_t drifts) override

Protected Member Functions inherited from img::DataProviderAlg
std::vector< float >	downscaleMax (std::size_t dst_size, std::vector< float > const &adc, size_t tick0) const

std::vector< float >	downscaleMaxMean (std::size_t dst_size, std::vector< float > const &adc, size_t tick0) const

std::vector< float >	downscaleMean (std::size_t dst_size, std::vector< float > const &adc, size_t tick0) const

std::vector< float >	downscale (std::size_t dst_size, std::vector< float > const &adc, size_t tick0) const

size_t	getDriftIndex (float drift) const

std::optional< std::vector< float > >	setWireData (std::vector< float > const &adc, size_t wireIdx) const

bool	patchFromDownsampledView (size_t wire, float drift, size_t size_w, size_t size_d, std::vector< std::vector< float >> &patch) const

bool	patchFromOriginalView (size_t wire, float drift, size_t size_w, size_t size_d, std::vector< std::vector< float >> &patch) const

Private Member Functions
WireDrift	getProjection (detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, const TLorentzVector &tvec, unsigned int plane) const

bool	setWireEdepsAndLabels (std::vector< float > const &edeps, std::vector< int > const &pdgs, size_t wireIdx)

void	collectVtxFlags (std::unordered_map< size_t, std::unordered_map< int, int >> &wireToDriftToVtxFlags, detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, const std::unordered_map< int, const simb::MCParticle * > &particleMap, unsigned int plane) const

bool	isElectronEnd (const simb::MCParticle &particle, const std::unordered_map< int, const simb::MCParticle * > &particleMap) const

bool	isMuonDecaying (const simb::MCParticle &particle, const std::unordered_map< int, const simb::MCParticle * > &particleMap) const

Static Private Member Functions
static float	particleRange2 (const simb::MCParticle &particle)

Private Attributes
double	fEdepTot

std::vector< std::vector< float > >	fWireDriftEdep

std::vector< std::vector< int > >	fWireDriftPdg

art::InputTag	fWireProducerLabel

art::InputTag	fHitProducerLabel

art::InputTag	fTrackModuleLabel

art::InputTag	fSimulationProducerLabel

art::InputTag	fSimChannelProducerLabel

bool	fSaveVtxFlags

bool	fSaveSimInfo

unsigned int	fAdcDelay

std::vector< size_t >	fEventsPerBin

Additional Inherited Members
Protected Attributes inherited from img::DataProviderAlg
DataProviderAlgView	fAlgView

EDownscaleMode	fDownscaleMode

size_t	fDriftWindow

bool	fDownscaleFullView

float	fDriftWindowInv

calo::CalorimetryAlg	fCalorimetryAlg

geo::GeometryCore const *	fGeometry

Detailed Description

Definition at line 198 of file PointIdAlg.h.

Member Enumeration Documentation

enum nnet::TrainingDataAlg::EMask

Enumerator
kNone
kPdgMask
kTypeMask
kVtxMask

Definition at line 200 of file PointIdAlg.h.

              {
     kNone = 0,
     kPdgMask = 0x00000FFF,  // pdg code mask
     kTypeMask = 0x0000F000, // track type mask
     kVtxMask = 0xFFFF0000   // vertex flags
   };

enum nnet::TrainingDataAlg::ETrkType

Enumerator
kDelta
kMichel
kPriEl
kPriMu

Definition at line 207 of file PointIdAlg.h.

                 {
     kDelta = 0x1000,  // delta electron
     kMichel = 0x2000, // Michel electron
     kPriEl = 0x4000,  // primary electron
     kPriMu = 0x8000   // primary muon
   };

enum nnet::TrainingDataAlg::EVtxId

Enumerator
kNuNC
kNuCC
kNuPri
kNuE
kNuMu
kNuTau
kHadr
kPi0
kDecay
kConv
kElectronEnd
kElastic
kInelastic

Definition at line 214 of file PointIdAlg.h.

               {
     kNuNC = 0x0010000,
     kNuCC = 0x0020000,
     kNuPri = 0x0040000, // nu interaction type
     kNuE = 0x0100000,
     kNuMu = 0x0200000,
     kNuTau = 0x0400000,        // nu flavor
     kHadr = 0x1000000,         // hadronic inelastic scattering
     kPi0 = 0x2000000,          // pi0 produced in this vertex
     kDecay = 0x4000000,        // point of particle decay
     kConv = 0x8000000,         // gamma conversion
     kElectronEnd = 0x10000000, // clear end of an electron
     kElastic = 0x20000000,     // Elastic scattering
     kInelastic = 0x40000000    // Inelastic scattering
   };

Constructor & Destructor Documentation

nnet::TrainingDataAlg::TrainingDataAlg ( const fhicl::ParameterSet & pset )

inline

Definition at line 254 of file PointIdAlg.h.

255 : TrainingDataAlg(fhicl::Table<Config>(pset, {})())

256 {}

fhicl::Table< Config >

nnet::TrainingDataAlg::TrainingDataAlg

TrainingDataAlg(const fhicl::ParameterSet &pset)

Definition: PointIdAlg.h:254

nnet::TrainingDataAlg::TrainingDataAlg ( const Config & config )

Definition at line 334 of file PointIdAlg.cxx.

   : img::DataProviderAlg(config)
   , fEdepTot(0)
   , fWireProducerLabel(config.WireLabel())
   , fHitProducerLabel(config.HitLabel())
   , fTrackModuleLabel(config.TrackLabel())
   , fSimulationProducerLabel(config.SimulationLabel())
   , fSimChannelProducerLabel(config.SimChannelLabel())
   , fSaveVtxFlags(config.SaveVtxFlags())
   , fAdcDelay(config.AdcDelayTicks())
   , fEventsPerBin(100, 0)
 {
   // If no sim channel producer is set then make it the same as the simulation label
   if(fSimChannelProducerLabel.label().empty())
     fSimChannelProducerLabel = fSimulationProducerLabel;
 
   fSaveSimInfo = !fSimulationProducerLabel.label().empty();
 }

nnet::TrainingDataAlg::~TrainingDataAlg ( )

overridedefault

Member Function Documentation

void nnet::TrainingDataAlg::collectVtxFlags	(	std::unordered_map< size_t, std::unordered_map< int, int >> &	wireToDriftToVtxFlags,
		detinfo::DetectorClocksData const &	clockData,
		detinfo::DetectorPropertiesData const &	detProp,
		const std::unordered_map< int, const simb::MCParticle * > &	particleMap,
		unsigned int	plane
	)		const

private

Definition at line 554 of file PointIdAlg.cxx.

 {
   for (auto const& p : particleMap) {
     auto const& particle = *p.second;
 
     double ekStart = 1000. * (particle.E() - particle.Mass());
     double ekEnd = 1000. * (particle.EndE() - particle.Mass());
     int pdg = abs(particle.PdgCode());
     int flagsStart = nnet::TrainingDataAlg::kNone;
     int flagsEnd = nnet::TrainingDataAlg::kNone;
 
     switch (pdg) {
     case 22:                                                     // gamma
       if ((particle.EndProcess() == "conv") && (ekStart > 40.0)) // conversion, gamma > 40MeV
       {
         flagsEnd = nnet::TrainingDataAlg::kConv;
       }
       break;
 
     case 11: // e+/-
       if (isElectronEnd(particle, particleMap)) { flagsEnd = nnet::TrainingDataAlg::kElectronEnd; }
       break;
 
     case 13: // mu+/-
       if (isMuonDecaying(particle, particleMap)) { flagsEnd = nnet::TrainingDataAlg::kDecay; }
       break;
 
     case 111: // pi0
       flagsStart = nnet::TrainingDataAlg::kPi0;
       break;
 
     case 321:  // K+/-
     case 211:  // pi+/-
     case 2212: // proton
       if (ekStart > 50.0) {
         if (particle.Mother() != 0) {
           auto search = particleMap.find(particle.Mother());
           if (search != particleMap.end()) {
             auto const& mother = *((*search).second);
             int m_pdg = abs(mother.PdgCode());
             unsigned int nSec = mother.NumberDaughters();
             unsigned int nVisible = 0;
             if (nSec > 1) {
               for (size_t d = 0; d < nSec; ++d) {
                 auto d_search = particleMap.find(mother.Daughter(d));
                 if (d_search != particleMap.end()) {
                   auto const& daughter = *((*d_search).second);
                   int d_pdg = abs(daughter.PdgCode());
                   if (((d_pdg == 2212) || (d_pdg == 211) || (d_pdg == 321)) &&
                       (1000. * (daughter.E() - daughter.Mass()) > 50.0)) {
                     ++nVisible;
                   }
                 }
               }
             }
             // hadron with Ek > 50MeV (so well visible) and
             // produced by another hadron (but not neutron, so not single track from nothing) or
             // at least secondary hadrons with Ek > 50MeV (so this is a good kink or V-like)
             if (((m_pdg != pdg) && (m_pdg != 2112)) || ((m_pdg != 2112) && (nVisible > 0)) ||
                 ((m_pdg == 2112) && (nVisible > 1))) {
               flagsStart = nnet::TrainingDataAlg::kHadr;
             }
           }
         }
 
         if (particle.EndProcess() == "FastScintillation") // potential decay at rest
         {
           unsigned int nSec = particle.NumberDaughters();
           for (size_t d = 0; d < nSec; ++d) {
             auto d_search = particleMap.find(particle.Daughter(d));
             if (d_search != particleMap.end()) {
               auto const& daughter = *((*d_search).second);
               int d_pdg = abs(daughter.PdgCode());
               if ((pdg == 321) && (d_pdg == 13)) {
                 flagsEnd = nnet::TrainingDataAlg::kDecay;
                 break;
               }
               if ((pdg == 211) && (d_pdg == 13)) {
                 flagsEnd = nnet::TrainingDataAlg::kDecay;
                 break;
               }
             }
           }
         }
 
         if ((particle.EndProcess() == "Decay") && (ekEnd > 200.0)) // decay in flight
         {
           unsigned int nSec = particle.NumberDaughters();
           for (size_t d = 0; d < nSec; ++d) {
             auto d_search = particleMap.find(particle.Daughter(d));
             if (d_search != particleMap.end()) {
               auto const& daughter = *((*d_search).second);
               int d_pdg = abs(daughter.PdgCode());
               if ((pdg == 321) && (d_pdg == 13)) {
                 flagsEnd = nnet::TrainingDataAlg::kHadr;
                 break;
               }
               if ((pdg == 211) && (d_pdg == 13)) {
                 flagsEnd = nnet::TrainingDataAlg::kHadr;
                 break;
               }
             }
           }
         }
       }
       break;
 
     default: continue;
     }
 
     if (particle.Process() == "primary") { flagsStart |= nnet::TrainingDataAlg::kNuPri; }
 
     if (flagsStart != nnet::TrainingDataAlg::kNone) {
       auto wd = getProjection(clockData, detProp, particle.Position(), plane);
 
       if ((wd.TPC == TPC()) && (wd.Cryo == Cryo())) {
         wireToDriftToVtxFlags[wd.Wire][wd.Drift] |= flagsStart;
       }
     }
     if (flagsEnd != nnet::TrainingDataAlg::kNone) {
       auto wd = getProjection(clockData, detProp, particle.EndPosition(), plane);
       if ((wd.TPC == TPC()) && (wd.Cryo == Cryo())) {
         wireToDriftToVtxFlags[wd.Wire][wd.Drift] |= flagsEnd;
       }
     }
 
     //TY: check elastic/inelastic scattering
     if (pdg == 321 || pdg == 211 || pdg == 2212) {
       simb::MCTrajectory truetraj = particle.Trajectory();
       auto thisTrajectoryProcessMap1 = truetraj.TrajectoryProcesses();
       if (thisTrajectoryProcessMap1.size()) {
         for (auto const& couple1 : thisTrajectoryProcessMap1) {
           if ((truetraj.KeyToProcess(couple1.second)).find("Elastic") != std::string::npos) {
             auto wd = getProjection(clockData, detProp, truetraj.at(couple1.first).first, plane);
             if ((wd.TPC == TPC()) && (wd.Cryo == Cryo())) {
               wireToDriftToVtxFlags[wd.Wire][wd.Drift] |= nnet::TrainingDataAlg::kElastic;
             }
           }
           if ((truetraj.KeyToProcess(couple1.second)).find("Inelastic") != std::string::npos) {
             auto wd = getProjection(clockData, detProp, truetraj.at(couple1.first).first, plane);
             if ((wd.TPC == TPC()) && (wd.Cryo == Cryo())) {
               wireToDriftToVtxFlags[wd.Wire][wd.Drift] |= nnet::TrainingDataAlg::kInelastic;
             }
           }
         }
       }
     }
   }
 }

bool nnet::TrainingDataAlg::findCrop	(	float	max_e_cut,
		unsigned int &	w0,
		unsigned int &	w1,
		unsigned int &	d0,
		unsigned int &	d1
	)		const

Definition at line 1020 of file PointIdAlg.cxx.

 {
   if (fWireDriftEdep.empty() || fWireDriftEdep.front().empty()) return false;
 
   float max_cut = 0.25 * max_e_cut;
 
   w0 = 0;
   float cut = 0;
   while (w0 < fWireDriftEdep.size()) {
     for (auto const d : fWireDriftEdep[w0])
       cut += d;
     if (cut < max_cut)
       w0++;
     else
       break;
   }
   w1 = fWireDriftEdep.size() - 1;
   cut = 0;
   while (w1 > w0) {
     for (auto const d : fWireDriftEdep[w1])
       cut += d;
     if (cut < max_cut)
       w1--;
     else
       break;
   }
   w1++;
 
   d0 = 0;
   cut = 0;
   while (d0 < fWireDriftEdep.front().size()) {
     for (size_t i = w0; i < w1; ++i)
       cut += fWireDriftEdep[i][d0];
     if (cut < max_cut)
       d0++;
     else
       break;
   }
   d1 = fWireDriftEdep.front().size() - 1;
   cut = 0;
   while (d1 > d0) {
     for (size_t i = w0; i < w1; ++i)
       cut += fWireDriftEdep[i][d1];
     if (cut < max_cut)
       d1--;
     else
       break;
   }
   d1++;
 
   unsigned int margin = 20;
   if ((w1 - w0 > 8) && (d1 - d0 > 8)) {
     if (w0 < margin)
       w0 = 0;
     else
       w0 -= margin;
 
     if (w1 > fWireDriftEdep.size() - margin)
       w1 = fWireDriftEdep.size();
     else
       w1 += margin;
 
     if (d0 < margin)
       d0 = 0;
     else
       d0 -= margin;
 
     if (d1 > fWireDriftEdep.front().size() - margin)
       d1 = fWireDriftEdep.front().size();
     else
       d1 += margin;
 
     return true;
   }
   else
     return false;
 }

double nnet::TrainingDataAlg::getEdepTot ( ) const

inline

Definition at line 293 of file PointIdAlg.h.

   {
     return fEdepTot;
   } // [GeV]

nnet::TrainingDataAlg::WireDrift nnet::TrainingDataAlg::getProjection	(	detinfo::DetectorClocksData const &	clockData,
		detinfo::DetectorPropertiesData const &	detProp,
		const TLorentzVector &	tvec,
		unsigned int	plane
	)		const

private

Definition at line 425 of file PointIdAlg.cxx.

 {
   nnet::TrainingDataAlg::WireDrift wd;
   wd.Wire = 0;
   wd.Drift = 0;
   wd.TPC = -1;
   wd.Cryo = -1;
 
   try {
     double vtx[3] = {tvec.X(), tvec.Y(), tvec.Z()};
     if (fGeometry->FindTPCAtPosition(vtx).isValid) {
       geo::TPCID tpcid = fGeometry->FindTPCAtPosition(vtx);
       unsigned int tpc = tpcid.TPC, cryo = tpcid.Cryostat;
 
       // correct for the time offset
       float dx = tvec.T() * 1.e-3 * detProp.DriftVelocity();
       int driftDir = fGeometry->TPC(tpcid).DetectDriftDirection();
       if (driftDir == 1) { dx *= -1; }
       else if (driftDir != -1) {
         throw cet::exception("nnet::TrainingDataAlg") << "drift direction is not X." << std::endl;
       }
       vtx[0] = tvec.X() + dx;
 
       wd.Wire = fGeometry->NearestWire(vtx, plane, tpc, cryo);
       wd.Drift = detProp.ConvertXToTicks(vtx[0], plane, tpc, cryo);
       wd.TPC = tpc;
       wd.Cryo = cryo;
     }
   }
   catch (const geo::InvalidWireIDError& e) {
     mf::LogWarning("TrainingDataAlg")
       << "Vertex projection out of wire planes, just skipping this vertex.";
   }
   catch (...) {
     mf::LogWarning("TrainingDataAlg") << "Vertex projection out of wire planes, skip MC vertex.";
   }
   return wd;
 }

bool nnet::TrainingDataAlg::isElectronEnd	(	const simb::MCParticle &	particle,
		const std::unordered_map< int, const simb::MCParticle * > &	particleMap
	)		const

private

Definition at line 469 of file PointIdAlg.cxx.

 {
   const float minElectronLength2 = 2.5 * 2.5;
   const float maxDeltaLength2 = 0.15 * 0.15;
 
   int pdg = abs(particle.PdgCode());
   if (pdg != 11) return false; // should be applied only to electrons
 
   size_t nSec = particle.NumberDaughters();
   for (size_t d = 0; d < nSec; ++d) {
     auto d_search = particleMap.find(particle.Daughter(d));
     if (d_search != particleMap.end()) {
       auto const& daughter = *((*d_search).second);
       int d_pdg = abs(daughter.PdgCode());
       if (d_pdg != 22) { return false; } // not the end of the shower
     }
   }
 
   float trkLength2 = 0;
   auto const* p = &particle;
   bool branching = false;
   while (!branching) {
     trkLength2 += particleRange2(*p);
     auto m_search = particleMap.find(p->Mother());
     if (m_search != particleMap.end()) {
       p = (*m_search).second;
       int m_pdg = abs(p->PdgCode());
       if (m_pdg == 11) {
         nSec = p->NumberDaughters();
         size_t ne = 0;
         for (size_t d = 0; d < nSec; ++d) {
           auto d_search = particleMap.find(p->Daughter(d));
           if (d_search != particleMap.end()) {
             auto const& daughter = *((*d_search).second);
             int d_pdg = abs(daughter.PdgCode());
             if (d_pdg == 11) {
               if (particleRange2(daughter) > maxDeltaLength2) { ne++; }
             }
           }
         }
         if (ne > 1) { branching = true; }
       }
       else
         break;
     }
     else
       break;
   }
 
   return (trkLength2 > minElectronLength2);
 }

bool nnet::TrainingDataAlg::isMuonDecaying	(	const simb::MCParticle &	particle,
		const std::unordered_map< int, const simb::MCParticle * > &	particleMap
	)		const

private

Definition at line 524 of file PointIdAlg.cxx.

 {
   bool hasElectron = false, hasNuMu = false, hasNuE = false;
 
   int pdg = abs(particle.PdgCode());
   //if ((pdg == 13) && (particle.EndProcess() == "FastScintillation" || particle.EndProcess() == "Decay" || particle.EndProcess() == "muMinusCaptureAtRest")) // potential muon decay at rest
   if ((pdg == 13) && (particle.EndProcess() == "FastScintillation" || particle.EndProcess() == "Decay")) // potential muon decay at rest
   {
     unsigned int nSec = particle.NumberDaughters();
     for (size_t d = 0; d < nSec; ++d) {
       auto d_search = particleMap.find(particle.Daughter(d));
       if (d_search != particleMap.end()) {
         auto const& daughter = *((*d_search).second);
         int d_pdg = abs(daughter.PdgCode());
         if (d_pdg == 11)
           hasElectron = true;
         else if (d_pdg == 14)
           hasNuMu = true;
         else if (d_pdg == 12)
           hasNuE = true;
       }
     }
   }
   
   return (hasElectron && hasNuMu && hasNuE);
 }

static float nnet::TrainingDataAlg::particleRange2 ( const simb::MCParticle & particle )

inlinestaticprivate

Definition at line 340 of file PointIdAlg.h.

   {
     float dx = particle.EndX() - particle.Vx();
     float dy = particle.EndY() - particle.Vy();
     float dz = particle.EndZ() - particle.Vz();
     return dx * dx + dy * dy + dz * dz;
   }

void nnet::TrainingDataAlg::reconfigure ( const Config & config )

img::DataProviderAlgView nnet::TrainingDataAlg::resizeView	(	detinfo::DetectorClocksData const &	clock_data,
		detinfo::DetectorPropertiesData const &	det_prop,
		size_t	wires,
		size_t	drifts
	)

overrideprotectedvirtual

Reimplemented from img::DataProviderAlg.

Definition at line 358 of file PointIdAlg.cxx.

 {
   auto view = img::DataProviderAlg::resizeView(clock_data, det_prop, wires, drifts);
 
   fWireDriftEdep.resize(wires);
   for (auto& w : fWireDriftEdep) {
     w.resize(view.fNCachedDrifts);
     std::fill(w.begin(), w.end(), 0.0F);
   }
 
   fWireDriftPdg.resize(wires);
   for (auto& w : fWireDriftPdg) {
     w.resize(view.fNCachedDrifts);
     std::fill(w.begin(), w.end(), 0);
   }
   return view;
 }

bool nnet::TrainingDataAlg::saveSimInfo ( ) const

inline

Definition at line 265 of file PointIdAlg.h.

   {
     return fSaveSimInfo;
   }

bool nnet::TrainingDataAlg::setDataEventData	(	const art::Event &	event,
		detinfo::DetectorClocksData const &	clockData,
		detinfo::DetectorPropertiesData const &	detProp,
		unsigned int	plane,
		unsigned int	tpc,
		unsigned int	cryo
	)

Definition at line 711 of file PointIdAlg.cxx.

 {
 
   art::Handle<std::vector<recob::Wire>> wireHandle;
   std::vector<art::Ptr<recob::Wire>> Wirelist;
 
   if (event.getByLabel(fWireProducerLabel, wireHandle)) art::fill_ptr_vector(Wirelist, wireHandle);
 
   if (!setWireDriftData(clockData, detProp, *wireHandle, plane, tpc, cryo)) {
     mf::LogError("TrainingDataAlg") << "Wire data not set.";
     return false;
   }
 
   // Hit info
   art::Handle<std::vector<recob::Hit>> HitHandle;
   std::vector<art::Ptr<recob::Hit>> Hitlist;
 
   if (event.getByLabel(fHitProducerLabel, HitHandle)) art::fill_ptr_vector(Hitlist, HitHandle);
 
   // Track info
   art::Handle<std::vector<recob::Track>> TrackHandle;
   std::vector<art::Ptr<recob::Track>> Tracklist;
 
   if (event.getByLabel(fTrackModuleLabel, TrackHandle))
     art::fill_ptr_vector(Tracklist, TrackHandle);
 
   art::FindManyP<recob::Track> ass_trk_hits(HitHandle, event, fTrackModuleLabel);
 
   // Loop over wires (sorry about hard coded value) to fill in 1) pdg and 2) charge depo
   for (size_t widx = 0; widx < 240; ++widx) {
 
     std::vector<float> labels_deposit(fAlgView.fNDrifts, 0); // full-drift-length buffers
     std::vector<int> labels_pdg(fAlgView.fNDrifts, 0);
 
     // First, the charge depo
     for (size_t subwidx = 0; subwidx < Wirelist.size(); ++subwidx) {
       if (widx + 240 == Wirelist[subwidx]->Channel()) {
         labels_deposit = Wirelist[subwidx]->Signal();
         break;
       }
     }
 
     // Second, the pdg code
     // This code finds the angle of the track and records
     //  events based on its angle to try to get an isometric sample
     //  instead of just a bunch of straight tracks
 
     // Meta code:
     // For each hit:
     //  find farthest hit from point
     //  then find farthest hit from THAT one
     //  should be start and end of track, then just use trig
 
     for (size_t iHit = 0; iHit < Hitlist.size(); ++iHit) {
 
       if (Hitlist[iHit]->Channel() != widx + 240) { continue; }
       if (Hitlist[iHit]->View() != 1) { continue; }
 
       // Make sure there is a track association
       if (ass_trk_hits.at(iHit).size() == 0) { continue; }
 
       // Not sure about this
       // Cutting on length to not just get a bunch of shower stubs
       // Might add a lot of bias though
       if (ass_trk_hits.at(iHit)[0]->Length() < 5) { continue; }
 
       // Search for farest hit from this one
       int far_index = 0;
       double far_dist = 0;
 
       for (size_t jHit = 0; jHit < Hitlist.size(); ++jHit) {
         if (jHit == iHit) { continue; }
         if (Hitlist[jHit]->View() != 1) { continue; }
 
         if (ass_trk_hits.at(jHit).size() == 0) { continue; }
         if (ass_trk_hits.at(jHit)[0]->ID() != ass_trk_hits.at(iHit)[0]->ID()) { continue; }
 
         double dist = sqrt((Hitlist[iHit]->Channel() - Hitlist[jHit]->Channel()) *
                              (Hitlist[iHit]->Channel() - Hitlist[jHit]->Channel()) +
                            (Hitlist[iHit]->PeakTime() - Hitlist[jHit]->PeakTime()) *
                              (Hitlist[iHit]->PeakTime() - Hitlist[jHit]->PeakTime()));
 
         if (far_dist < dist) {
           far_dist = dist;
           far_index = jHit;
         }
       }
 
       // Search for the other end of the track
       int other_end = 0;
       int other_dist = 0;
 
       for (size_t jHit = 0; jHit < Hitlist.size(); ++jHit) {
         if (jHit == iHit or int(jHit) == far_index) { continue; }
         if (Hitlist[jHit]->View() != 1) { continue; }
 
         if (ass_trk_hits.at(jHit).size() == 0) { continue; }
         if (ass_trk_hits.at(jHit)[0]->ID() != ass_trk_hits.at(iHit)[0]->ID()) { continue; }
 
         double dist = sqrt((Hitlist[far_index]->Channel() - Hitlist[jHit]->Channel()) *
                              (Hitlist[far_index]->Channel() - Hitlist[jHit]->Channel()) +
                            (Hitlist[far_index]->PeakTime() - Hitlist[jHit]->PeakTime()) *
                              (Hitlist[far_index]->PeakTime() - Hitlist[jHit]->PeakTime()));
 
         if (other_dist < dist) {
           other_dist = dist;
           other_end = jHit;
         }
       }
 
       // We have the end points now
       double del_wire = double(Hitlist[other_end]->Channel() - Hitlist[far_index]->Channel());
       double del_time = double(Hitlist[other_end]->PeakTime() - Hitlist[far_index]->PeakTime());
       double hypo = sqrt(del_wire * del_wire + del_time * del_time);
 
       if (hypo == 0) { continue; } // Should never happen, but doing it anyway
 
       double cosser = TMath::Abs(del_wire / hypo);
       double norm_ang = TMath::ACos(cosser) * 2 / TMath::Pi();
 
       // Using fEventsPerBin to keep track of number of hits per angle (normalized to 0 to 1)
 
       int binner = int(norm_ang * fEventsPerBin.size());
       if (binner >= (int)fEventsPerBin.size()) {
         binner = fEventsPerBin.size() - 1;
       } // Dealing with rounding errors
 
       // So we should get a total of 5000 * 100 = 50,000 if we use the whole set
       if (fEventsPerBin[binner] > 5000) { continue; }
       fEventsPerBin[binner]++;
 
       // If survives everything, saves the pdg
       labels_pdg[Hitlist[iHit]->PeakTime()] = 211; // Same as pion for now
     }
 
     setWireEdepsAndLabels(labels_deposit, labels_pdg, widx);
 
   } // for each Wire
 
   return true;
 }

bool nnet::TrainingDataAlg::setEventData	(	const art::Event &	event,
		detinfo::DetectorClocksData const &	clockData,
		detinfo::DetectorPropertiesData const &	detProp,
		unsigned int	plane,
		unsigned int	tpc,
		unsigned int	cryo
	)

Definition at line 859 of file PointIdAlg.cxx.

 {
   art::ValidHandle<std::vector<recob::Wire>> wireHandle =
     event.getValidHandle<std::vector<recob::Wire>>(fWireProducerLabel);
 
   if (!setWireDriftData(clockData, detProp, *wireHandle, plane, tpc, cryo)) {
     mf::LogError("TrainingDataAlg") << "Wire data not set.";
     return false;
   }
 
   if (!fSaveSimInfo || event.isRealData()) {
     mf::LogInfo("TrainingDataAlg") << "Skip MC simulation info.";
     return true;
   }
 
   art::ServiceHandle<sim::LArG4Parameters const> larParameters;
   double electronsToGeV = 1. / larParameters->GeVToElectrons();
 
   auto particleHandle =
     event.getValidHandle<std::vector<simb::MCParticle>>(fSimulationProducerLabel);
 
   auto simChannelHandle =
     event.getValidHandle<std::vector<sim::SimChannel>>(fSimChannelProducerLabel);
 
   std::unordered_map<int, const simb::MCParticle*> particleMap;
   for (auto const& particle : *particleHandle) {
     particleMap[particle.TrackId()] = &particle;
   }
 
   std::unordered_map<size_t, std::unordered_map<int, int>> wireToDriftToVtxFlags;
   if (fSaveVtxFlags) collectVtxFlags(wireToDriftToVtxFlags, clockData, detProp, particleMap, plane);
 
   fEdepTot = 0;
 
   std::map<int, int> trackToPDG;
   for (size_t widx = 0; widx < fAlgView.fNWires; ++widx) {
     auto wireChannelNumber = fAlgView.fWireChannels[widx];
     if (wireChannelNumber == raw::InvalidChannelID) continue;
 
     std::vector<float> labels_deposit(fAlgView.fNDrifts, 0); // full-drift-length buffers,
     std::vector<int> labels_pdg(labels_deposit.size(), 0);   // both of the same size,
     int labels_size = labels_deposit.size();                 // cached as int for comparisons below
 
     std::map<int, std::map<int, double>> timeToTrackToCharge;
     for (auto const& channel : *simChannelHandle) {
       if (channel.Channel() != wireChannelNumber) continue;
 
       auto const& timeSlices = channel.TDCIDEMap();
       for (auto const& timeSlice : timeSlices) {
         int time = timeSlice.first;
 
         auto const& energyDeposits = timeSlice.second;
         for (auto const& energyDeposit : energyDeposits) {
           int pdg = 0;
           int tid = energyDeposit.trackID;
           if (tid < 0) // negative tid means it is EM activity, and -tid is the mother
           {
             pdg = 11;
             tid = -tid;
 
             auto search = particleMap.find(tid);
             if (search == particleMap.end()) {
               mf::LogWarning("TrainingDataAlg") << "PARTICLE NOT FOUND";
               continue;
             }
             auto const& mother = *((*search).second); // mother particle of this EM
             int mPdg = abs(mother.PdgCode());
             if ((mPdg == 13) || (mPdg == 211) || (mPdg == 2212)) {
               if (energyDeposit.numElectrons > 10)
                 pdg |= nnet::TrainingDataAlg::kDelta; // tag delta ray
             }
           }
           else {
             auto search = particleMap.find(tid);
             if (search == particleMap.end()) {
               mf::LogWarning("TrainingDataAlg") << "PARTICLE NOT FOUND";
               continue;
             }
             auto const& particle = *((*search).second);
             pdg = abs(particle.PdgCode());
 
             if (particle.Process() == "primary") {
               if (pdg == 11) {
                 pdg |= nnet::TrainingDataAlg::kPriEl; // tag primary
               }
               else if (pdg == 13) {
                 pdg |= nnet::TrainingDataAlg::kPriMu; // tag primary
               }
             }
 
             /*
             auto msearch = particleMap.find(particle.Mother());
             if (msearch != particleMap.end()) {
               auto const& mother = *((*msearch).second);
               if (pdg == 11) // electron, check if it is Michel
               {
                 if (nnet::TrainingDataAlg::isMuonDecaying(mother, particleMap)) {
                   std::cout<<particle.Process()<<std::endl;
                   pdg |= nnet::TrainingDataAlg::kMichel; // tag Michel
                 }
               }
             }
             */
             if (pdg == 11){ // electron, check if it is Michel or delta ray
               if (particle.Process() == "Decay"){
                 pdg |= nnet::TrainingDataAlg::kMichel; // tag Michel
               }
               else if (particle.Process() == "muIoni"){
                 pdg |= nnet::TrainingDataAlg::kDelta; // tag delta ray
               }
             }
           }
 
           trackToPDG[energyDeposit.trackID] = pdg;
 
           double energy = energyDeposit.numElectrons * electronsToGeV;
           timeToTrackToCharge[time][energyDeposit.trackID] += energy;
           fEdepTot += energy;
 
         } // loop over energy deposits
       }   // loop over time slices
     }     // for each SimChannel
 
     int type_pdg_mask = nnet::TrainingDataAlg::kTypeMask | nnet::TrainingDataAlg::kPdgMask;
     for (auto const& ttc : timeToTrackToCharge) {
       float max_deposit = 0.0;
       int max_pdg = 0;
       for (auto const& tc : ttc.second) {
 
         if (tc.second > max_deposit) {
           max_deposit = tc.second;
           max_pdg = trackToPDG[tc.first];
         }
       }
 
       if (ttc.first < labels_size) {
         int tick_idx = ttc.first + fAdcDelay;
         if (tick_idx < labels_size) {
           labels_deposit[tick_idx] = max_deposit;
           labels_pdg[tick_idx] = max_pdg & type_pdg_mask;
         }
       }
     }
 
     for (auto const& drift_flags : wireToDriftToVtxFlags[widx]) {
       int drift = drift_flags.first, flags = drift_flags.second;
       if ((drift >= 0) && (drift < labels_size)) { labels_pdg[drift] |= flags; }
     }
     setWireEdepsAndLabels(labels_deposit, labels_pdg, widx);
   } // for each Wire
 
   return true;
 }

bool nnet::TrainingDataAlg::setWireEdepsAndLabels	(	std::vector< float > const &	edeps,
		std::vector< int > const &	pdgs,
		size_t	wireIdx
	)

private

Definition at line 381 of file PointIdAlg.cxx.

 {
   if ((wireIdx >= fWireDriftEdep.size()) || (edeps.size() != pdgs.size())) { return false; }
 
   size_t dstep = 1;
   if (fDownscaleFullView) { dstep = fDriftWindow; }
 
   if (edeps.size() / dstep > fAlgView.fNCachedDrifts) { return false; }
 
   auto& wEdep = fWireDriftEdep[wireIdx];
   auto& wPdg = fWireDriftPdg[wireIdx];
 
   for (size_t i = 0; i < fAlgView.fNCachedDrifts; ++i) {
     size_t i0 = i * dstep;
     size_t i1 = (i + 1) * dstep;
 
     int best_pdg = pdgs[i0] & nnet::TrainingDataAlg::kPdgMask;
     int vtx_flags = pdgs[i0] & nnet::TrainingDataAlg::kVtxMask;
     int type_flags = pdgs[i0] & nnet::TrainingDataAlg::kTypeMask;
     float max_edep = edeps[i0];
     for (size_t k = i0 + 1; k < i1; ++k) {
       float ek = edeps[k];
       if (ek > max_edep) {
         max_edep = ek;
         best_pdg = pdgs[k] & nnet::TrainingDataAlg::kPdgMask; // remember best matching pdg
       }
       type_flags |= pdgs[k] & nnet::TrainingDataAlg::kTypeMask; // accumulate track type flags
       vtx_flags |= pdgs[k] & nnet::TrainingDataAlg::kVtxMask;   // accumulate all vtx flags
     }
 
     wEdep[i] = max_edep;
 
     best_pdg |= type_flags;
     if (fSaveVtxFlags) best_pdg |= vtx_flags;
     wPdg[i] = best_pdg;
   }
 
   return true;
 }

std::vector<float> const& nnet::TrainingDataAlg::wireEdep ( size_t widx ) const

inline

Definition at line 298 of file PointIdAlg.h.

   {
     return fWireDriftEdep[widx];
   }

std::vector<int> const& nnet::TrainingDataAlg::wirePdg ( size_t widx ) const

inline

Definition at line 303 of file PointIdAlg.h.

   {
     return fWireDriftPdg[widx];
   }

Member Data Documentation

unsigned int nnet::TrainingDataAlg::fAdcDelay

private

Definition at line 365 of file PointIdAlg.h.

double nnet::TrainingDataAlg::fEdepTot

private

Definition at line 353 of file PointIdAlg.h.

std::vector<size_t> nnet::TrainingDataAlg::fEventsPerBin

private

Definition at line 367 of file PointIdAlg.h.

art::InputTag nnet::TrainingDataAlg::fHitProducerLabel

private

Definition at line 358 of file PointIdAlg.h.

bool nnet::TrainingDataAlg::fSaveSimInfo

private

Definition at line 363 of file PointIdAlg.h.

bool nnet::TrainingDataAlg::fSaveVtxFlags

private

Definition at line 362 of file PointIdAlg.h.

art::InputTag nnet::TrainingDataAlg::fSimChannelProducerLabel

private

Definition at line 361 of file PointIdAlg.h.

art::InputTag nnet::TrainingDataAlg::fSimulationProducerLabel

private

Definition at line 360 of file PointIdAlg.h.

art::InputTag nnet::TrainingDataAlg::fTrackModuleLabel

private

Definition at line 359 of file PointIdAlg.h.

std::vector<std::vector<float> > nnet::TrainingDataAlg::fWireDriftEdep

private

Definition at line 354 of file PointIdAlg.h.

std::vector<std::vector<int> > nnet::TrainingDataAlg::fWireDriftPdg

private

Definition at line 355 of file PointIdAlg.h.

art::InputTag nnet::TrainingDataAlg::fWireProducerLabel

private

Definition at line 357 of file PointIdAlg.h.

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

larrecodnn/larrecodnn/ImagePatternAlgs/Tensorflow/PointIdAlg/PointIdAlg.h
larrecodnn/larrecodnn/ImagePatternAlgs/Tensorflow/PointIdAlg/PointIdAlg.cxx

Classes

Public Types

Public Member Functions

Protected Member Functions

Private Member Functions

Static Private Member Functions

Private Attributes

Additional Inherited Members

Detailed Description

Member Enumeration Documentation

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation