Classes
class	AverageWaveform

class	CalibrationAnalysis

class	DefaultOpDetResponse

class	DUNE35tonOpDetResponse

class	DUNEOpDetResponse

class	DUNEOpDetResponseInterface

class	FIFOHistogramAna

class	FlashHypothesis

class	FlashHypothesisAna

class	FlashHypothesisAnaAlg

class	FlashHypothesisCalculator

class	FlashHypothesisCollection

class	FlashHypothesisComparison

class	FlashHypothesisCreator

class	FlashMatchAna

class	FlashUtilities

class	FocusList

class	LEDCalibrationAna

class	MicrobooneOpDetResponse

class	OpDetDigiAnaDUNE

class	OpDetDigitizerDUNE

class	OpDetDigitizerDUNEDP

class	OpDetDigitizerProtoDUNE

class	OpDetResponseInterface

class	OpDigiAna

class	OpDigiProperties

class	OpFlashAna

class	OpFlashAnaAlg

class	OpFlashFinder

class	OpFlashFinderDualPhase

class	OpFlashMCTruthAna

class	OpFlashSimpleAna

class	OpHitAna

class	OpHitFinder

class	OpMCDigi

class	OpSlicer

class	OptDetDigitizer

class	OpticalRawDigitReformatter

class	PDSNoiseFilter

class	ProtoDUNE_opdet_eventdisplay

class	ProtoDUNEOpDetResponse

class	SimPhotonCounter

class	SimPhotonCounterAlg

class	SIPMOpSensorSim

class	WaveformDigitizerSim

Typedefs
typedef std::vector< std::pair< size_t, size_t > >	Ranges_t

Functions
bool	sortOpHitByTime (const art::Ptr< recob::OpHit > &left, const art::Ptr< recob::OpHit > &right)

void	writeHistogram (std::vector< double > const &binned)

void	checkOnBeamFlash (std::vector< recob::OpFlash > const &FlashVector)

void	RunFlashFinder (std::vector< recob::OpHit > const &HitVector, std::vector< recob::OpFlash > &FlashVector, std::vector< std::vector< int >> &AssocList, double const BinWidth, geo::GeometryCore const &geom, float const FlashThreshold, float const WidthTolerance, detinfo::DetectorClocksData const &ClocksData, float const TrigCoinc)

unsigned int	GetAccumIndex (double const PeakTime, double const MinTime, double const BinWidth, double const BinOffset)

void	FillAccumulator (unsigned int const &AccumIndex, unsigned int const &HitIndex, double const PE, float const FlashThreshold, std::vector< double > &Binned, std::vector< std::vector< int >> &Contributors, std::vector< int > &FlashesInAccumulator)

void	FillFlashesBySizeMap (std::vector< int > const &FlashesInAccumulator, std::vector< double > const &BinnedPE, int const &Accumulator, std::map< double, std::map< int, std::vector< int >>, std::greater< double >> &FlashesBySize)

void	FillHitsThisFlash (std::vector< std::vector< int >> const &Contributors, int const &Bin, std::vector< int > const &HitClaimedByFlash, std::vector< int > &HitsThisFlash)

void	ClaimHits (std::vector< recob::OpHit > const &HitVector, std::vector< int > const &HitsThisFlash, float const FlashThreshold, std::vector< std::vector< int >> &HitsPerFlash, std::vector< int > &HitClaimedByFlash)

void	AssignHitsToFlash (std::vector< int > const &FlashesInAccumulator1, std::vector< int > const &FlashesInAccumulator2, std::vector< double > const &Binned1, std::vector< double > const &Binned2, std::vector< std::vector< int >> const &Contributors1, std::vector< std::vector< int >> const &Contributors2, std::vector< recob::OpHit > const &HitVector, std::vector< std::vector< int >> &HitsPerFlash, float const FlashThreshold)

void	FindSeedHit (std::map< double, std::vector< int >, std::greater< double >> const &HitsBySize, std::vector< bool > &HitsUsed, std::vector< recob::OpHit > const &HitVector, std::vector< int > &HitsThisRefinedFlash, double &PEAccumulated, double &FlashMaxTime, double &FlashMinTime)

void	AddHitToFlash (int const &HitID, std::vector< bool > &HitsUsed, recob::OpHit const &currentHit, double const WidthTolerance, std::vector< int > &HitsThisRefinedFlash, double &PEAccumulated, double &FlashMaxTime, double &FlashMinTime)

void	CheckAndStoreFlash (std::vector< std::vector< int >> &RefinedHitsPerFlash, std::vector< int > const &HitsThisRefinedFlash, double const PEAccumulated, float const FlashThreshold, std::vector< bool > &HitsUsed)

void	RefineHitsInFlash (std::vector< int > const &HitsThisFlash, std::vector< recob::OpHit > const &HitVector, std::vector< std::vector< int >> &RefinedHitsPerFlash, float const WidthTolerance, float const FlashThreshold)

void	AddHitContribution (recob::OpHit const &currentHit, double &MaxTime, double &MinTime, double &AveTime, double &FastToTotal, double &AveAbsTime, double &TotalPE, std::vector< double > &PEs)

void	GetHitGeometryInfo (recob::OpHit const &currentHit, geo::GeometryCore const &geom, std::vector< double > &sumw, std::vector< double > &sumw2, double &sumy, double &sumy2, double &sumz, double &sumz2)

double	CalculateWidth (double const sum, double const sum_squared, double const weights_sum)

void	ConstructFlash (std::vector< int > const &HitsPerFlashVec, std::vector< recob::OpHit > const &HitVector, std::vector< recob::OpFlash > &FlashVector, geo::GeometryCore const &geom, detinfo::DetectorClocksData const &ClocksData, float const TrigCoinc)

double	GetLikelihoodLateLight (double const iPE, double const iTime, double const iWidth, double const jPE, double const jTime, double const jWidth)

void	MarkFlashesForRemoval (std::vector< recob::OpFlash > const &FlashVector, size_t const BeginFlash, std::vector< bool > &MarkedForRemoval)

void	RemoveFlashesFromVectors (std::vector< bool > const &MarkedForRemoval, std::vector< recob::OpFlash > &FlashVector, size_t const BeginFlash, std::vector< std::vector< int >> &RefinedHitsPerFlash)

void	RemoveLateLight (std::vector< recob::OpFlash > &FlashVector, std::vector< std::vector< int >> &RefinedHitsPerFlash)

template<typename T , typename Compare >
std::vector< int >	sort_permutation (std::vector< T > const &vec, int offset, Compare compare)

template<typename T >
void	apply_permutation (std::vector< T > &vec, std::vector< int > const &p)

void	RunHitFinder (std::vector< raw::OpDetWaveform > const &opDetWaveformVector, std::vector< recob::OpHit > &hitVector, pmtana::PulseRecoManager const &pulseRecoMgr, pmtana::PMTPulseRecoBase const &threshAlg, geo::GeometryCore const &geometry, float hitThreshold, detinfo::DetectorClocksData const &clocksData, calib::IPhotonCalibrator const &calibrator, bool use_start_time)

void	ConstructHit (float hitThreshold, int channel, double timeStamp, pmtana::pulse_param const &pulse, std::vector< recob::OpHit > &hitVector, detinfo::DetectorClocksData const &clocksData, calib::IPhotonCalibrator const &calibrator, bool use_start_time)

Detailed Description

Title: FlashHypothesis Class Author: Wes Ketchum (wketc.nosp@m.hum@.nosp@m.lanl..nosp@m.gov)

Description: Class that represents a flash hypothesis (PEs per opdet) and its error

Title: FlashHypothesis Calculator Class Author: Wes Ketchum (wketc.nosp@m.hum@.nosp@m.lanl..nosp@m.gov)

Description: Simple class for calculating flash hypotheses

Title: FlashHypothesis Creator Class Author: Wes Ketchum (wketc.nosp@m.hum@.nosp@m.lanl..nosp@m.gov)

Description: class that produces a flash hypothesis for a trajectory. Input: Trajectory (std::vector<TVector3> objects) Output: FlashHypotheses

Title: FlashUtilities Class Author: Wes Ketchum (wketc.nosp@m.hum@.nosp@m.lanl..nosp@m.gov)

Description: Class that contains utility functions for flash and flash hypotheses: — compare a flash hypothesis to a truth or reco vector — get an extent of a flash (central point, width) These classes should operate using simple objects, and will need other classes/functions to fill those vectors properly.

Title: OpFlash Algorithims Authors, editors: Ben Jones, MIT Wes Ketchum wketc.nosp@m.hum@.nosp@m.lanl..nosp@m.gov Gleb Sinev gleb..nosp@m.sine.nosp@m.v@duk.nosp@m.e.ed.nosp@m.u Alex Himmel ahimm.nosp@m.el@f.nosp@m.nal.g.nosp@m.ov

Description: These are the algorithms used by OpFlashFinder to produce flashes.

Title: OpHit Algorithims Authors, editors: Ben Jones, MIT Wes Ketchum wketc.nosp@m.hum@.nosp@m.lanl..nosp@m.gov Gleb Sinev gleb..nosp@m.sine.nosp@m.v@duk.nosp@m.e.ed.nosp@m.u Alex Himmel ahimm.nosp@m.el@f.nosp@m.nal.g.nosp@m.ov Kevin Wood kevin.nosp@m..woo.nosp@m.d@sto.nosp@m.nybr.nosp@m.ook.e.nosp@m.du

Description: These are the algorithms used by OpHitFinder to produce optical hits.

Title: SimPhotonCounter Class Author: Wes Ketchum (wketc.nosp@m.hum@.nosp@m.lanl..nosp@m.gov)

Description: Class that counts up sim photons, leading towards comparisons with flashes and flash hypotheses.

Title: SimPhotonCounterALG Class Author: Wes Ketchum (wketc.nosp@m.hum@.nosp@m.lanl..nosp@m.gov)

Description: Alg class that counts up sim photons, leading towards comparisons with flashes and flash hypotheses.

Typedef Documentation

typedef std::vector< std::pair< size_t, size_t > > opdet::Ranges_t

Definition at line 75 of file WaveformDigitizerSim_module.cc.

Function Documentation

void opdet::AddHitContribution	(	recob::OpHit const &	currentHit,
		double &	MaxTime,
		double &	MinTime,
		double &	AveTime,
		double &	FastToTotal,
		double &	AveAbsTime,
		double &	TotalPE,
		std::vector< double > &	PEs
	)

Definition at line 455 of file OpFlashAlg.cxx.

   {
     double PEThisHit = currentHit.PE();
     double TimeThisHit = currentHit.PeakTime();
     if (TimeThisHit > MaxTime) MaxTime = TimeThisHit;
     if (TimeThisHit < MinTime) MinTime = TimeThisHit;
 
     // These quantities for the flash are defined
     // as the weighted averages over the hits
     AveTime += PEThisHit * TimeThisHit;
     FastToTotal += PEThisHit * currentHit.FastToTotal();
     AveAbsTime += PEThisHit * currentHit.PeakTimeAbs();
 
     // These are totals
     TotalPE += PEThisHit;
     PEs.at(static_cast<unsigned int>(currentHit.OpChannel())) += PEThisHit;
   }

void opdet::AddHitToFlash	(	int const &	HitID,
		std::vector< bool > &	HitsUsed,
		recob::OpHit const &	currentHit,
		double const	WidthTolerance,
		std::vector< int > &	HitsThisRefinedFlash,
		double &	PEAccumulated,
		double &	FlashMaxTime,
		double &	FlashMinTime
	)

Definition at line 330 of file OpFlashAlg.cxx.

   {
     if (HitsUsed.at(HitID)) return;
 
     double HitTime = currentHit.PeakTime();
     double HitWidth = 0.5 * currentHit.Width();
     double FlashTime = 0.5 * (FlashMaxTime + FlashMinTime);
     double FlashWidth = 0.5 * (FlashMaxTime - FlashMinTime);
 
     if (std::abs(HitTime - FlashTime) > WidthTolerance * (HitWidth + FlashWidth)) return;
 
     HitsThisRefinedFlash.push_back(HitID);
     FlashMaxTime = std::max(FlashMaxTime, HitTime + HitWidth);
     FlashMinTime = std::min(FlashMinTime, HitTime - HitWidth);
     PEAccumulated += currentHit.PE();
     HitsUsed[HitID] = true;
 
   } // End AddHitToFlash

template<typename T >

void opdet::apply_permutation	(	std::vector< T > &	vec,
		std::vector< int > const &	p
	)

Definition at line 703 of file OpFlashAlg.cxx.

   {
 
     std::vector<T> sorted_vec(p.size());
     std::transform(p.begin(), p.end(), sorted_vec.begin(), [&](int i) { return vec[i]; });
     vec = sorted_vec;
   }

void opdet::AssignHitsToFlash	(	std::vector< int > const &	FlashesInAccumulator1,
		std::vector< int > const &	FlashesInAccumulator2,
		std::vector< double > const &	Binned1,
		std::vector< double > const &	Binned2,
		std::vector< std::vector< int >> const &	Contributors1,
		std::vector< std::vector< int >> const &	Contributors2,
		std::vector< recob::OpHit > const &	HitVector,
		std::vector< std::vector< int >> &	HitsPerFlash,
		float const	FlashThreshold
	)

Definition at line 249 of file OpFlashAlg.cxx.

   {
     // Sort all the flashes found by size. The structure is:
     // FlashesBySize[flash size][accumulator_num] = [flash_index1, flash_index2...]
     std::map<double, std::map<int, std::vector<int>>, std::greater<double>> FlashesBySize;
 
     // Sort the flashes by size using map
     FillFlashesBySizeMap(FlashesInAccumulator1, Binned1, 1, FlashesBySize);
     FillFlashesBySizeMap(FlashesInAccumulator2, Binned2, 2, FlashesBySize);
 
     // This keeps track of which hits are claimed by which flash
     std::vector<int> HitClaimedByFlash(HitVector.size(), -1);
 
     // Walk from largest to smallest, claiming hits.
     // The biggest flash always gets dibbs,
     // but we keep track of overlaps for re-merging later (do we? ---WK)
     for (auto const& itFlash : FlashesBySize)
       // If several with same size, walk through accumulators
       for (auto const& itAcc : itFlash.second) {
 
         int Accumulator = itAcc.first;
 
         // Walk through flash-tagged bins in this accumulator
         for (auto const& Bin : itAcc.second) {
 
           std::vector<int> HitsThisFlash;
 
           if (Accumulator == 1)
             FillHitsThisFlash(Contributors1, Bin, HitClaimedByFlash, HitsThisFlash);
           else if (Accumulator == 2)
             FillHitsThisFlash(Contributors2, Bin, HitClaimedByFlash, HitsThisFlash);
 
           ClaimHits(HitVector, HitsThisFlash, FlashThreshold, HitsPerFlash, HitClaimedByFlash);
 
         } // End loop over this accumulator
 
       } // End loops over accumulators
     // End of loops over sorted flashes
 
   } // End AssignHitsToFlash

double opdet::CalculateWidth	(	double const	sum,
		double const	sum_squared,
		double const	weights_sum
	)

Definition at line 514 of file OpFlashAlg.cxx.

   {
     if (sum_squared * weights_sum - sum * sum < 0)
       return 0;
     else
       return std::sqrt(sum_squared * weights_sum - sum * sum) / weights_sum;
   }

void opdet::CheckAndStoreFlash	(	std::vector< std::vector< int >> &	RefinedHitsPerFlash,
		std::vector< int > const &	HitsThisRefinedFlash,
		double const	PEAccumulated,
		float const	FlashThreshold,
		std::vector< bool > &	HitsUsed
	)

Definition at line 358 of file OpFlashAlg.cxx.

   {
     // If above threshold, we just add hits to the flash vector, and move on
     if (PEAccumulated >= FlashThreshold) {
       RefinedHitsPerFlash.push_back(HitsThisRefinedFlash);
       return;
     }
 
     // If there is only one hit in collection, we can immediately move on
     if (HitsThisRefinedFlash.size() == 1) return;
 
     // We need to release all other hits (allow possible reuse)
     for (std::vector<int>::const_iterator hitIterator = std::next(HitsThisRefinedFlash.begin());
          hitIterator != HitsThisRefinedFlash.end();
          ++hitIterator)
       HitsUsed.at(*hitIterator) = false;
 
   } // End CheckAndStoreFlash

void opdet::checkOnBeamFlash ( std::vector< recob::OpFlash > const & FlashVector )

Definition at line 53 of file OpFlashAlg.cxx.

   {
     for (auto const& flash : FlashVector)
       if (flash.OnBeamTime() == 1)
         std::cout << "OnBeamFlash with time " << flash.Time() << std::endl;
   }

void opdet::ClaimHits	(	std::vector< recob::OpHit > const &	HitVector,
		std::vector< int > const &	HitsThisFlash,
		float const	FlashThreshold,
		std::vector< std::vector< int >> &	HitsPerFlash,
		std::vector< int > &	HitClaimedByFlash
	)

Definition at line 226 of file OpFlashAlg.cxx.

   {
     // Check for newly claimed hits
     double PE = 0;
     for (auto const& Hit : HitsThisFlash)
       PE += HitVector.at(Hit).PE();
 
     if (PE < FlashThreshold) return;
 
     // Add the flash to the list
     HitsPerFlash.push_back(HitsThisFlash);
 
     // And claim all the hits
     for (auto const& Hit : HitsThisFlash)
       if (HitClaimedByFlash.at(Hit) == -1) HitClaimedByFlash.at(Hit) = HitsPerFlash.size() - 1;
   }

void opdet::ConstructFlash	(	std::vector< int > const &	HitsPerFlashVec,
		std::vector< recob::OpHit > const &	HitVector,
		std::vector< recob::OpFlash > &	FlashVector,
		geo::GeometryCore const &	geom,
		detinfo::DetectorClocksData const &	ClocksData,
		float const	TrigCoinc
	)

Definition at line 524 of file OpFlashAlg.cxx.

   {
     double MaxTime = -std::numeric_limits<double>::max();
     double MinTime = std::numeric_limits<double>::max();
 
     std::vector<double> PEs(geom.MaxOpChannel() + 1, 0.0);
     unsigned int Nplanes = geom.Nplanes();
     std::vector<double> sumw(Nplanes, 0.0);
     std::vector<double> sumw2(Nplanes, 0.0);
 
     double TotalPE = 0;
     double AveTime = 0;
     double AveAbsTime = 0;
     double FastToTotal = 0;
     double sumy = 0;
     double sumz = 0;
     double sumy2 = 0;
     double sumz2 = 0;
 
     for (auto const& HitID : HitsPerFlashVec) {
       AddHitContribution(
         HitVector.at(HitID), MaxTime, MinTime, AveTime, FastToTotal, AveAbsTime, TotalPE, PEs);
       GetHitGeometryInfo(HitVector.at(HitID), geom, sumw, sumw2, sumy, sumy2, sumz, sumz2);
     }
 
     AveTime /= TotalPE;
     AveAbsTime /= TotalPE;
     FastToTotal /= TotalPE;
 
     double meany = sumy / TotalPE;
     double meanz = sumz / TotalPE;
 
     double widthy = CalculateWidth(sumy, sumy2, TotalPE);
     double widthz = CalculateWidth(sumz, sumz2, TotalPE);
 
     std::vector<double> WireCenters(Nplanes, 0.0);
     std::vector<double> WireWidths(Nplanes, 0.0);
 
     for (size_t p = 0; p != Nplanes; ++p) {
       WireCenters.at(p) = sumw.at(p) / TotalPE;
       WireWidths.at(p) = CalculateWidth(sumw.at(p), sumw2.at(p), TotalPE);
     }
 
     // Emprical corrections to get the Frame right.
     // Eventual solution - remove frames
     int Frame = ClocksData.OpticalClock().Frame(AveAbsTime - 18.1);
     if (Frame == 0) Frame = 1;
 
     int BeamFrame = ClocksData.OpticalClock().Frame(ClocksData.TriggerTime());
     bool InBeamFrame = false;
     if (!(ClocksData.TriggerTime() < 0)) InBeamFrame = (Frame == BeamFrame);
 
     double TimeWidth = (MaxTime - MinTime) / 2.0;
 
     int OnBeamTime = 0;
     if (InBeamFrame && (std::abs(AveTime) < TrigCoinc)) OnBeamTime = 1;
 
     FlashVector.emplace_back(AveTime,
                              TimeWidth,
                              AveAbsTime,
                              Frame,
                              PEs,
                              InBeamFrame,
                              OnBeamTime,
                              FastToTotal,
                              meany,
                              widthy,
                              meanz,
                              widthz,
                              WireCenters,
                              WireWidths);
   }

void opdet::ConstructHit	(	float	hitThreshold,
		int	channel,
		double	timeStamp,
		pmtana::pulse_param const &	pulse,
		std::vector< recob::OpHit > &	hitVector,
		detinfo::DetectorClocksData const &	clocksData,
		calib::IPhotonCalibrator const &	calibrator,
		bool	use_start_time
	)

Definition at line 64 of file OpHitAlg.cxx.

   {
 
     if (pulse.peak < hitThreshold) return;
 
     double absTime = timeStamp + clocksData.OpticalClock().TickPeriod() * (use_start_time ? pulse.t_start : pulse.t_max);
 
     double relTime = absTime - clocksData.TriggerTime();
 
     int frame = clocksData.OpticalClock().Frame(timeStamp);
 
     double PE = 0.0;
     if (calibrator.UseArea())
       PE = calibrator.PE(pulse.area, channel);
     else
       PE = calibrator.PE(pulse.peak, channel);
 
     double width = (pulse.t_end - pulse.t_start) * clocksData.OpticalClock().TickPeriod();
 
     hitVector.emplace_back(
       channel, relTime, absTime, frame, width, pulse.area, pulse.peak, PE, 0.0);
   }

void opdet::FillAccumulator	(	unsigned int const &	AccumIndex,
		unsigned int const &	HitIndex,
		double const	PE,
		float const	FlashThreshold,
		std::vector< double > &	Binned,
		std::vector< std::vector< int >> &	Contributors,
		std::vector< int > &	FlashesInAccumulator
	)

Definition at line 181 of file OpFlashAlg.cxx.

   {
 
     Contributors.at(AccumIndex).push_back(HitIndex);
 
     Binned.at(AccumIndex) += PE;
 
     // If this wasn't a flash already, add it to the list
     if (Binned.at(AccumIndex) >= FlashThreshold && (Binned.at(AccumIndex) - PE) < FlashThreshold)
       FlashesInAccumulator.push_back(AccumIndex);
   }

void opdet::FillFlashesBySizeMap	(	std::vector< int > const &	FlashesInAccumulator,
		std::vector< double > const &	BinnedPE,
		int const &	Accumulator,
		std::map< double, std::map< int, std::vector< int >>, std::greater< double >> &	FlashesBySize
	)

Definition at line 201 of file OpFlashAlg.cxx.

   {
     for (auto const& flash : FlashesInAccumulator)
       FlashesBySize[BinnedPE.at(flash)][Accumulator].push_back(flash);
   }

void opdet::FillHitsThisFlash	(	std::vector< std::vector< int >> const &	Contributors,
		int const &	Bin,
		std::vector< int > const &	HitClaimedByFlash,
		std::vector< int > &	HitsThisFlash
	)

Definition at line 213 of file OpFlashAlg.cxx.

   {
     // For each hit in the flash
     for (auto const& HitIndex : Contributors.at(Bin))
       // If unclaimed, claim it
       if (HitClaimedByFlash.at(HitIndex) == -1) HitsThisFlash.push_back(HitIndex);
   }

void opdet::FindSeedHit	(	std::map< double, std::vector< int >, std::greater< double >> const &	HitsBySize,
		std::vector< bool > &	HitsUsed,
		std::vector< recob::OpHit > const &	HitVector,
		std::vector< int > &	HitsThisRefinedFlash,
		double &	PEAccumulated,
		double &	FlashMaxTime,
		double &	FlashMinTime
	)

Definition at line 300 of file OpFlashAlg.cxx.

   {
     for (auto const& itHit : HitsBySize)
       for (auto const& HitID : itHit.second) {
 
         if (HitsUsed.at(HitID)) continue;
 
         PEAccumulated = HitVector.at(HitID).PE();
         FlashMaxTime = HitVector.at(HitID).PeakTime() + 0.5 * HitVector.at(HitID).Width();
         FlashMinTime = HitVector.at(HitID).PeakTime() - 0.5 * HitVector.at(HitID).Width();
 
         HitsThisRefinedFlash.clear();
         HitsThisRefinedFlash.push_back(HitID);
 
         HitsUsed.at(HitID) = true;
         return;
 
       } // End loop over inner vector
     // End loop over HitsBySize map
 
   } // End FindSeedHit

unsigned int opdet::GetAccumIndex	(	double const	PeakTime,
		double const	MinTime,
		double const	BinWidth,
		double const	BinOffset
	)

Definition at line 171 of file OpFlashAlg.cxx.

   {
     return static_cast<unsigned int>((PeakTime - MinTime + BinOffset) / BinWidth);
   }

void opdet::GetHitGeometryInfo	(	recob::OpHit const &	currentHit,
		geo::GeometryCore const &	geom,
		std::vector< double > &	sumw,
		std::vector< double > &	sumw2,
		double &	sumy,
		double &	sumy2,
		double &	sumz,
		double &	sumz2
	)

Definition at line 482 of file OpFlashAlg.cxx.

   {
     double xyz[3];
     geom.OpDetGeoFromOpChannel(currentHit.OpChannel()).GetCenter(xyz);
     double PEThisHit = currentHit.PE();
 
     geo::TPCID tpc = geom.FindTPCAtPosition(xyz);
     // if the point does not fall into any TPC,
     // it does not contribute to the average wire position
     if (tpc.isValid) {
       for (size_t p = 0; p != geom.Nplanes(); ++p) {
         geo::PlaneID const planeID(tpc, p);
         unsigned int w = geom.NearestWire(xyz, planeID);
         sumw.at(p) += PEThisHit * w;
         sumw2.at(p) += PEThisHit * w * w;
       }
     } // if we found the TPC
     sumy += PEThisHit * xyz[1];
     sumy2 += PEThisHit * xyz[1] * xyz[1];
     sumz += PEThisHit * xyz[2];
     sumz2 += PEThisHit * xyz[2] * xyz[2];
   }

double opdet::GetLikelihoodLateLight	(	double const	iPE,
		double const	iTime,
		double const	iWidth,
		double const	jPE,
		double const	jTime,
		double const	jWidth
	)

Definition at line 604 of file OpFlashAlg.cxx.

   {
     if (iTime > jTime) return 1e6;
 
     // Calculate hypothetical PE if this were actually a late flash from i.
     // Argon time const is 1600 ns, so 1.6.
     double HypPE = iPE * jWidth / iWidth * std::exp(-(jTime - iTime) / 1.6);
     double nsigma = (jPE - HypPE) / std::sqrt(HypPE);
     return nsigma;
   }

void opdet::MarkFlashesForRemoval	(	std::vector< recob::OpFlash > const &	FlashVector,
		size_t const	BeginFlash,
		std::vector< bool > &	MarkedForRemoval
	)

Definition at line 622 of file OpFlashAlg.cxx.

   {
     for (size_t iFlash = BeginFlash; iFlash != FlashVector.size(); ++iFlash) {
 
       double iTime = FlashVector.at(iFlash).Time();
       double iPE = FlashVector.at(iFlash).TotalPE();
       double iWidth = FlashVector.at(iFlash).TimeWidth();
 
       for (size_t jFlash = iFlash + 1; jFlash != FlashVector.size(); ++jFlash) {
 
         if (MarkedForRemoval.at(jFlash - BeginFlash)) continue;
 
         double jTime = FlashVector.at(jFlash).Time();
         double jPE = FlashVector.at(jFlash).TotalPE();
         double jWidth = FlashVector.at(jFlash).TimeWidth();
 
         // If smaller than, or within 2sigma of expectation,
         // attribute to late light and toss out
         if (GetLikelihoodLateLight(iPE, iTime, iWidth, jPE, jTime, jWidth) < 3.0)
           MarkedForRemoval.at(jFlash - BeginFlash) = true;
       }
     }
   }

void opdet::RefineHitsInFlash	(	std::vector< int > const &	HitsThisFlash,
		std::vector< recob::OpHit > const &	HitVector,
		std::vector< std::vector< int >> &	RefinedHitsPerFlash,
		float const	WidthTolerance,
		float const	FlashThreshold
	)

Definition at line 383 of file OpFlashAlg.cxx.

   {
     // Sort the hits by their size using map
     // HitsBySize[HitSize] = [hit1, hit2 ...]
     std::map<double, std::vector<int>, std::greater<double>> HitsBySize;
     for (auto const& HitID : HitsThisFlash)
       HitsBySize[HitVector.at(HitID).PE()].push_back(HitID);
 
     // Heres what we do:
     //  1.Start with the biggest remaining hit
     //  2.Look for any within one width of this hit
     //  3.Find the new upper and lower bounds of the flash
     //  4.Collect again
     //  5.Repeat until no new hits collected
     //  6.Remove these hits from consideration and repeat
 
     std::vector<bool> HitsUsed(HitVector.size(), false);
     double PEAccumulated, FlashMaxTime, FlashMinTime;
     std::vector<int> HitsThisRefinedFlash;
 
     while (true) {
 
       HitsThisRefinedFlash.clear();
       PEAccumulated = 0;
       FlashMaxTime = 0;
       FlashMinTime = 0;
 
       FindSeedHit(HitsBySize,
                   HitsUsed,
                   HitVector,
                   HitsThisRefinedFlash,
                   PEAccumulated,
                   FlashMaxTime,
                   FlashMinTime);
 
       if (HitsThisRefinedFlash.size() == 0) return;
 
       // Start this at zero to do the while at least once
       size_t NHitsThisRefinedFlash = 0;
 
       // If size of HitsThisRefinedFlash is same size,
       // that means we're not adding anymore
       while (NHitsThisRefinedFlash < HitsThisRefinedFlash.size()) {
         NHitsThisRefinedFlash = HitsThisRefinedFlash.size();
 
         for (auto const& itHit : HitsBySize)
           for (auto const& HitID : itHit.second)
             AddHitToFlash(HitID,
                           HitsUsed,
                           HitVector.at(HitID),
                           WidthTolerance,
                           HitsThisRefinedFlash,
                           PEAccumulated,
                           FlashMaxTime,
                           FlashMinTime);
       }
 
       // We did our collecting, now check if the flash is
       // still good and push back
       CheckAndStoreFlash(
         RefinedHitsPerFlash, HitsThisRefinedFlash, PEAccumulated, FlashThreshold, HitsUsed);
 
     } // End while there are hits left
 
   } // End RefineHitsInFlash

void opdet::RemoveFlashesFromVectors	(	std::vector< bool > const &	MarkedForRemoval,
		std::vector< recob::OpFlash > &	FlashVector,
		size_t const	BeginFlash,
		std::vector< std::vector< int >> &	RefinedHitsPerFlash
	)

Definition at line 650 of file OpFlashAlg.cxx.

   {
     for (int iFlash = MarkedForRemoval.size() - 1; iFlash != -1; --iFlash)
       if (MarkedForRemoval.at(iFlash)) {
         RefinedHitsPerFlash.erase(RefinedHitsPerFlash.begin() + iFlash);
         FlashVector.erase(FlashVector.begin() + BeginFlash + iFlash);
       }
   }

void opdet::RemoveLateLight	(	std::vector< recob::OpFlash > &	FlashVector,
		std::vector< std::vector< int >> &	RefinedHitsPerFlash
	)

Definition at line 664 of file OpFlashAlg.cxx.

   {
     std::vector<bool> MarkedForRemoval(RefinedHitsPerFlash.size(), false);
 
     size_t const BeginFlash = FlashVector.size() - RefinedHitsPerFlash.size();
 
     recob::OpFlashSortByTime sort_flash_by_time;
 
     // Determine the sort of FlashVector starting at BeginFlash
     auto sort_order = sort_permutation(FlashVector, BeginFlash, sort_flash_by_time);
 
     // Sort the RefinedHitsPerFlash in the same way as tail end of FlashVector
     apply_permutation(RefinedHitsPerFlash, sort_order);
 
     std::sort(FlashVector.begin() + BeginFlash, FlashVector.end(), sort_flash_by_time);
 
     MarkFlashesForRemoval(FlashVector, BeginFlash, MarkedForRemoval);
 
     RemoveFlashesFromVectors(MarkedForRemoval, FlashVector, BeginFlash, RefinedHitsPerFlash);
 
   } // End RemoveLateLight

void opdet::RunFlashFinder	(	std::vector< recob::OpHit > const &	HitVector,
		std::vector< recob::OpFlash > &	FlashVector,
		std::vector< std::vector< int >> &	AssocList,
		double const	BinWidth,
		geo::GeometryCore const &	geom,
		float const	FlashThreshold,
		float const	WidthTolerance,
		detinfo::DetectorClocksData const &	ClocksData,
		float const	TrigCoinc
	)

Definition at line 62 of file OpFlashAlg.cxx.

   {
     // Initial size for accumulators - will be automatically extended if needed
     int initialsize = 6400;
 
     // These are the accumulators which will hold broad-binned light yields
     std::vector<double> Binned1(initialsize);
     std::vector<double> Binned2(initialsize);
 
     // These will keep track of which pulses put activity in each bin
     std::vector<std::vector<int>> Contributors1(initialsize);
     std::vector<std::vector<int>> Contributors2(initialsize);
 
     // These will keep track of where we have met the flash condition
     // (in order to prevent second pointless loop)
     std::vector<int> FlashesInAccumulator1;
     std::vector<int> FlashesInAccumulator2;
 
     double minTime = std::numeric_limits<float>::max();
     for (auto const& hit : HitVector)
       if (hit.PeakTime() < minTime) minTime = hit.PeakTime();
 
     for (auto const& hit : HitVector) {
 
       double peakTime = hit.PeakTime();
 
       unsigned int AccumIndex1 = GetAccumIndex(peakTime, minTime, BinWidth, 0.0);
 
       unsigned int AccumIndex2 = GetAccumIndex(peakTime, minTime, BinWidth, BinWidth / 2.0);
 
       // Extend accumulators if needed (2 always larger than 1)
       if (AccumIndex2 >= Binned1.size()) {
         std::cout << "Extending vectors to " << AccumIndex2 * 1.2 << std::endl;
         Binned1.resize(AccumIndex2 * 1.2);
         Binned2.resize(AccumIndex2 * 1.2);
         Contributors1.resize(AccumIndex2 * 1.2);
         Contributors2.resize(AccumIndex2 * 1.2);
       }
 
       size_t const hitIndex = &hit - &HitVector[0];
 
       FillAccumulator(AccumIndex1,
                       hitIndex,
                       hit.PE(),
                       FlashThreshold,
                       Binned1,
                       Contributors1,
                       FlashesInAccumulator1);
 
       FillAccumulator(AccumIndex2,
                       hitIndex,
                       hit.PE(),
                       FlashThreshold,
                       Binned2,
                       Contributors2,
                       FlashesInAccumulator2);
 
     } // End loop over hits
 
     // Now start to create flashes.
     // First, need vector to keep track of which hits belong to which flashes
     std::vector<std::vector<int>> HitsPerFlash;
 
     //if (Frame == 1) writeHistogram(Binned1);
 
     AssignHitsToFlash(FlashesInAccumulator1,
                       FlashesInAccumulator2,
                       Binned1,
                       Binned2,
                       Contributors1,
                       Contributors2,
                       HitVector,
                       HitsPerFlash,
                       FlashThreshold);
 
     // Now we do the fine grained part.
     // Subdivide each flash into sub-flashes with overlaps within hit widths
     // (assumed wider than photon travel time)
     std::vector<std::vector<int>> RefinedHitsPerFlash;
     for (auto const& HitsThisFlash : HitsPerFlash)
       RefineHitsInFlash(
         HitsThisFlash, HitVector, RefinedHitsPerFlash, WidthTolerance, FlashThreshold);
 
     // Now we have all our hits assigned to a flash.
     // Make the recob::OpFlash objects
     for (auto const& HitsPerFlashVec : RefinedHitsPerFlash)
       ConstructFlash(HitsPerFlashVec, HitVector, FlashVector, geom, ClocksData, TrigCoinc);
 
     RemoveLateLight(FlashVector, RefinedHitsPerFlash);
 
     //checkOnBeamFlash(FlashVector);
 
     // Finally, write the association list.
     // back_inserter tacks the result onto the end of AssocList
     for (auto& HitIndicesThisFlash : RefinedHitsPerFlash)
       AssocList.push_back(HitIndicesThisFlash);
 
   } // End RunFlashFinder

void opdet::RunHitFinder	(	std::vector< raw::OpDetWaveform > const &	opDetWaveformVector,
		std::vector< recob::OpHit > &	hitVector,
		pmtana::PulseRecoManager const &	pulseRecoMgr,
		pmtana::PMTPulseRecoBase const &	threshAlg,
		geo::GeometryCore const &	geometry,
		float	hitThreshold,
		detinfo::DetectorClocksData const &	clocksData,
		calib::IPhotonCalibrator const &	calibrator,
		bool	use_start_time
	)

Definition at line 29 of file OpHitAlg.cxx.

   {
 
     for (auto const& waveform : opDetWaveformVector) {
 
       const int channel = static_cast<int>(waveform.ChannelNumber());
 
       if (!geometry.IsValidOpChannel(channel)) {
         mf::LogError("OpHitFinder")
           << "Error! unrecognized channel number " << channel << ". Ignoring pulse";
         continue;
       }
 
       pulseRecoMgr.Reconstruct(waveform);
 
       // Get the result
       auto const& pulses = threshAlg.GetPulses();
 
       const double timeStamp = waveform.TimeStamp();
 
       for (auto const& pulse : pulses)
         ConstructHit(hitThreshold, channel, timeStamp, pulse, hitVector, clocksData, calibrator, use_start_time);
     }
   }

template<typename T , typename Compare >

std::vector< int > opdet::sort_permutation	(	std::vector< T > const &	vec,
		int	offset,
		Compare	compare
	)

Definition at line 690 of file OpFlashAlg.cxx.

   {
 
     std::vector<int> p(vec.size() - offset);
     std::iota(p.begin(), p.end(), 0);
     std::sort(
       p.begin(), p.end(), [&](int i, int j) { return compare(vec[i + offset], vec[j + offset]); });
     return p;
   }

bool opdet::sortOpHitByTime	(	const art::Ptr< recob::OpHit > &	left,
		const art::Ptr< recob::OpHit > &	right
	)

Definition at line 44 of file OpSlicer_module.cc.

   {
     return left->PeakTime() < right->PeakTime();
   }

void opdet::writeHistogram ( std::vector< double > const & binned )

Definition at line 34 of file OpFlashAlg.cxx.

   {
     TH1D* binned_histogram = new TH1D("binned_histogram",
                                       "Collection of All OpHits;Time (ms);PEs",
                                       binned.size(),
                                       0,
                                       binned.size());
     for (size_t i = 0; i < binned.size(); ++i)
       binned_histogram->SetBinContent(i, binned.at(i));
 
     TFile f_out("output_hist.root", "RECREATE");
     binned_histogram->Write();
     f_out.Close();
 
     delete binned_histogram;
   }

Classes

Typedefs

Functions

Detailed Description

Typedef Documentation

Function Documentation