Inheritance diagram for SpacePointAnalysisMC::SpacePointAnalysisMC:

Classes
struct	ideCompare

Public Member Functions
	SpacePointAnalysisMC (fhicl::ParameterSet const &pset)
	Constructor. More...

	~SpacePointAnalysisMC ()
	Destructor. More...

void	configure (fhicl::ParameterSet const &pset) override

void	initializeHists (art::ServiceHandle< art::TFileService > &, const std::string &) override
	Interface for initializing the histograms to be filled. More...

void	initializeTuple (TTree *) override
	Interface for initializing the tuple variables. More...

void	endJob (int numEvents) override
	Interface for method to executve at the end of run processing. More...

void	fillHistograms (const art::Event &) const override
	Interface for filling histograms. More...

Public Member Functions inherited from IHitEfficiencyHistogramTool
virtual	~IHitEfficiencyHistogramTool () noexcept=default
	Virtual Destructor. More...

virtual	~IHitEfficiencyHistogramTool () noexcept=default
	Virtual Destructor. More...

Private Types
using	SimIDESet = std::set< const sim::IDE *, ideCompare >

using	IDEToVoxelIDMap = std::unordered_map< const sim::IDE *, sim::LArVoxelID >

using	VoxelIDToIDESetMap = std::map< sim::LArVoxelID, SimIDESet >

using	TDCToIDEMap = std::map< unsigned short, SimIDESet >

using	ChanToTDCToIDEMap = std::map< raw::ChannelID_t, TDCToIDEMap >

using	VoxelIDSet = std::set< sim::LArVoxelID >

using	TDCToIDESetMap = std::unordered_map< unsigned short, SimIDESet >

using	PlaneToTDCToIDESetMap = std::map< unsigned short, TDCToIDESetMap >

using	VoxelIDToPlaneTDCIDEMap = std::map< sim::LArVoxelID, PlaneToTDCToIDESetMap >

using	TDCIDEPair = std::pair< unsigned short, const sim::IDE * >

using	TickTDCIDEVec = std::vector< TDCIDEPair >

using	ChanToTDCIDEMap = std::unordered_map< raw::ChannelID_t, TickTDCIDEVec >

using	TrackIDChanToTDCIDEMap = std::unordered_map< int, ChanToTDCIDEMap >

using	ChargeDeposit = std::tuple< TDCIDEPair, TDCIDEPair, TDCIDEPair, float, float >

using	ChargeDepositVec = std::vector< ChargeDeposit >

using	ChanToChargeMap = std::map< raw::ChannelID_t, ChargeDepositVec >

using	TrackToChanChargeMap = std::unordered_map< int, ChanToChargeMap >

using	HitPointerVec = std::vector< const recob::Hit * >

using	RecobHitToVoxelIDMap = std::unordered_map< const recob::Hit *, VoxelIDSet >

using	HitSimObjVec = std::vector< HitSimulationTupleObj >

Private Member Functions
void	clear () const

void	makeTrackToChanChargeMap (const TrackIDChanToTDCIDEMap &, TrackToChanChargeMap &, float &, int &) const

void	compareHitsToSim (const art::Event &, const ChanToTDCToIDEMap &, const ChanToChargeMap &, const ChanToTDCIDEMap &, const IDEToVoxelIDMap &, RecobHitToVoxelIDMap &) const

void	matchHitSim (const detinfo::DetectorClocksData &clockData, const HitPointerVec &, const ChanToTDCToIDEMap &, const ChargeDepositVec &, const ChanToTDCIDEMap &, const IDEToVoxelIDMap &, RecobHitToVoxelIDMap &) const

void	compareSpacePointsToSim (const art::Event &, const detinfo::DetectorClocksData &clockData, const RecobHitToVoxelIDMap &, const VoxelIDToPlaneTDCIDEMap &) const

Private Attributes
std::vector< art::InputTag >	fRecobHitLabelVec

std::vector< art::InputTag >	fSpacePointLabelVec

art::InputTag	fMCParticleProducerLabel

art::InputTag	fSimChannelProducerLabel

art::InputTag	fSimEnergyProducerLabel

art::InputTag	fBadChannelProducerLabel

std::vector< int >	fOffsetVec
	Allow offsets for each plane. More...

float	fSimChannelMinEnergy

float	fSimEnergyMinEnergy

TTree *	fTree

std::vector< int >	fTPCVec

std::vector< int >	fCryoVec

std::vector< int >	fPlaneVec

HitSimObjVec	fHitSimObjVec

HitSpacePointObj	fHitSpacePointObj

const geo::GeometryCore *	fGeometry
	pointer to Geometry service More...

Detailed Description

Definition at line 311 of file SpacePointAnalysisMC_tool.cc.

Member Typedef Documentation

using SpacePointAnalysisMC::SpacePointAnalysisMC::ChanToChargeMap = std::map<raw::ChannelID_t,ChargeDepositVec>

private

Definition at line 389 of file SpacePointAnalysisMC_tool.cc.

using SpacePointAnalysisMC::SpacePointAnalysisMC::ChanToTDCIDEMap = std::unordered_map<raw::ChannelID_t,TickTDCIDEVec>

private

Definition at line 383 of file SpacePointAnalysisMC_tool.cc.

using SpacePointAnalysisMC::SpacePointAnalysisMC::ChanToTDCToIDEMap = std::map<raw::ChannelID_t, TDCToIDEMap>

private

Definition at line 372 of file SpacePointAnalysisMC_tool.cc.

using SpacePointAnalysisMC::SpacePointAnalysisMC::ChargeDeposit = std::tuple<TDCIDEPair,TDCIDEPair,TDCIDEPair,float,float>

private

Definition at line 387 of file SpacePointAnalysisMC_tool.cc.

using SpacePointAnalysisMC::SpacePointAnalysisMC::ChargeDepositVec = std::vector<ChargeDeposit>

private

Definition at line 388 of file SpacePointAnalysisMC_tool.cc.

using SpacePointAnalysisMC::SpacePointAnalysisMC::HitPointerVec = std::vector<const recob::Hit*>

private

Definition at line 396 of file SpacePointAnalysisMC_tool.cc.

using SpacePointAnalysisMC::SpacePointAnalysisMC::HitSimObjVec = std::vector<HitSimulationTupleObj>

private

Definition at line 426 of file SpacePointAnalysisMC_tool.cc.

using SpacePointAnalysisMC::SpacePointAnalysisMC::IDEToVoxelIDMap = std::unordered_map<const sim::IDE*, sim::LArVoxelID>

private

Definition at line 369 of file SpacePointAnalysisMC_tool.cc.

using SpacePointAnalysisMC::SpacePointAnalysisMC::PlaneToTDCToIDESetMap = std::map<unsigned short, TDCToIDESetMap>

private

Definition at line 377 of file SpacePointAnalysisMC_tool.cc.

using SpacePointAnalysisMC::SpacePointAnalysisMC::RecobHitToVoxelIDMap = std::unordered_map<const recob::Hit*, VoxelIDSet>

private

Definition at line 397 of file SpacePointAnalysisMC_tool.cc.

using SpacePointAnalysisMC::SpacePointAnalysisMC::SimIDESet = std::set<const sim::IDE*,ideCompare>

private

Definition at line 368 of file SpacePointAnalysisMC_tool.cc.

using SpacePointAnalysisMC::SpacePointAnalysisMC::TDCIDEPair = std::pair<unsigned short, const sim::IDE*>

private

Definition at line 381 of file SpacePointAnalysisMC_tool.cc.

using SpacePointAnalysisMC::SpacePointAnalysisMC::TDCToIDEMap = std::map<unsigned short, SimIDESet>

private

Definition at line 371 of file SpacePointAnalysisMC_tool.cc.

using SpacePointAnalysisMC::SpacePointAnalysisMC::TDCToIDESetMap = std::unordered_map<unsigned short,SimIDESet>

private

Definition at line 376 of file SpacePointAnalysisMC_tool.cc.

using SpacePointAnalysisMC::SpacePointAnalysisMC::TickTDCIDEVec = std::vector<TDCIDEPair>

private

Definition at line 382 of file SpacePointAnalysisMC_tool.cc.

using SpacePointAnalysisMC::SpacePointAnalysisMC::TrackIDChanToTDCIDEMap = std::unordered_map<int,ChanToTDCIDEMap>

private

Definition at line 384 of file SpacePointAnalysisMC_tool.cc.

using SpacePointAnalysisMC::SpacePointAnalysisMC::TrackToChanChargeMap = std::unordered_map<int,ChanToChargeMap>

private

Definition at line 390 of file SpacePointAnalysisMC_tool.cc.

using SpacePointAnalysisMC::SpacePointAnalysisMC::VoxelIDSet = std::set<sim::LArVoxelID>

private

Definition at line 373 of file SpacePointAnalysisMC_tool.cc.

using SpacePointAnalysisMC::SpacePointAnalysisMC::VoxelIDToIDESetMap = std::map<sim::LArVoxelID, SimIDESet>

private

Definition at line 370 of file SpacePointAnalysisMC_tool.cc.

using SpacePointAnalysisMC::SpacePointAnalysisMC::VoxelIDToPlaneTDCIDEMap = std::map<sim::LArVoxelID, PlaneToTDCToIDESetMap>

private

Definition at line 378 of file SpacePointAnalysisMC_tool.cc.

Constructor & Destructor Documentation

SpacePointAnalysisMC::SpacePointAnalysisMC::SpacePointAnalysisMC ( fhicl::ParameterSet const & pset )

explicit

Constructor.

Parameters

pset	Constructor.

Arguments:

pset - Fcl parameters.

Definition at line 442 of file SpacePointAnalysisMC_tool.cc.

                                                                          : fTree(nullptr)
 {
     fGeometry           = lar::providerFrom<geo::Geometry>();
 
     configure(pset);
 
     // Report.
     mf::LogInfo("SpacePointAnalysisMC") << "SpacePointAnalysisMC configured\n";
 }

SpacePointAnalysisMC::SpacePointAnalysisMC::~SpacePointAnalysisMC ( )

Destructor.

Definition at line 454 of file SpacePointAnalysisMC_tool.cc.

455 {}

Member Function Documentation

void SpacePointAnalysisMC::SpacePointAnalysisMC::clear ( ) const

private

Definition at line 520 of file SpacePointAnalysisMC_tool.cc.

 {
     fTPCVec.clear();
     fCryoVec.clear();
     fPlaneVec.clear();
 
     fHitSpacePointObj.clear();
 
     for(auto& hitObj : fHitSimObjVec) hitObj.clear();
 
     return;
 }

void SpacePointAnalysisMC::SpacePointAnalysisMC::compareHitsToSim	(	const art::Event &	event,
		const ChanToTDCToIDEMap &	chanToTDCToIDEMap,
		const ChanToChargeMap &	chanToChargeMap,
		const ChanToTDCIDEMap &	chanToTDCIDEMap,
		const IDEToVoxelIDMap &	ideToVoxelIDMap,
		RecobHitToVoxelIDMap &	recobHitToVoxelIDMap
	)		const

private

Definition at line 701 of file SpacePointAnalysisMC_tool.cc.

 {
     // We start by building a mapping between channels and lists of hits on that channel (ordered by time)
     using HitPointerVec   = std::vector<const recob::Hit*>;
     using ChanToHitVecMap = std::map<raw::ChannelID_t,HitPointerVec>;
 
     ChanToHitVecMap chanToHitVecMap;
 
     // And now fill it
     for(const auto& hitLabel : fRecobHitLabelVec)
     {
         art::Handle< std::vector<recob::Hit> > hitHandle;
         event.getByLabel(hitLabel, hitHandle);
 
         // If no hits then skip
         if ((*hitHandle).empty()) continue;
 
         for(const auto& hit : *hitHandle) chanToHitVecMap[hit.Channel()].emplace_back(&hit);
     }
 
     // Now go through and order each vector of hits by time
     for(auto& chanToHitPair : chanToHitVecMap)
     {
         HitPointerVec& hitPtrVec = chanToHitPair.second;
 
         std::sort(hitPtrVec.begin(),
                   hitPtrVec.end(),
                   [](const auto& left, const auto& right){return left->Channel() == right->Channel() ? left->PeakTime() < right->PeakTime() : left->Channel() < right->Channel();});
     }
 
     auto const clockData = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(event);
 
     // The idea is to loop over the input sim information so we can look at efficiency as well as resolution issues
     for(const auto& chanToChargePair : chanToChargeMap)
     {
         // Recover the channel
         raw::ChannelID_t channel = chanToChargePair.first;
 
         // Look up the hits associated to this channel
         ChanToHitVecMap::const_iterator chanToHitVecItr = chanToHitVecMap.find(channel);
 
         // For now we simply punt...
         if (chanToHitVecItr == chanToHitVecMap.end()) continue;
 
         // Recover channel information based on this hit
         const ChargeDepositVec& chargeDepositVec = chanToChargePair.second;
 
         // Get the hits...
         const HitPointerVec& hitPtrVec = chanToHitVecItr->second;
 
         short int lastSnippetStart(-1);
 
         HitPointerVec hitVec;
 
         // Outer loop over hits in this hit collection
         for(const auto& hitPtr : hitPtrVec)
         {
             // We want to collect together the hits that are on the same snippet. Hits will come grouped and in order
             // along the snippet, so we simply keep them in a local vector until we hit the end of the snippet...
             //
             // ** It is worth noting this scheme as implemented will miss the last snippet of hits... so think about
             // that for the future
             short int snippetStart = hitPtr->StartTick();
 
             if (snippetStart == lastSnippetStart)
             {
                 lastSnippetStart = snippetStart;
 
                 hitVec.emplace_back(hitPtr);
 
                 continue;
             }
 
             // Process the current list of hits (which will be on the same snippet)
             matchHitSim(clockData, hitVec, chanToTDCToIDEMap, chargeDepositVec, chanToTDCIDEMap, ideToVoxelIDMap, recobHitToVoxelIDMap);
 
             hitVec.clear();
             hitVec.emplace_back(hitPtr);
             lastSnippetStart = snippetStart;
         }
 
         // Make sure to catch the last set of hits in the group
         if (!hitVec.empty()) matchHitSim(clockData, hitVec, chanToTDCToIDEMap, chargeDepositVec, chanToTDCIDEMap, ideToVoxelIDMap, recobHitToVoxelIDMap);
     }
 
     return;
 }

void SpacePointAnalysisMC::SpacePointAnalysisMC::compareSpacePointsToSim	(	const art::Event &	event,
		const detinfo::DetectorClocksData &	clockData,
		const RecobHitToVoxelIDMap &	recobHitToVoxelIDMap,
		const VoxelIDToPlaneTDCIDEMap &	voxelToPlaneTDCIDEMap
	)		const

private

Definition at line 960 of file SpacePointAnalysisMC_tool.cc.

 {
     // Armed with these maps we can now process the SpacePoints...
     if (!recobHitToVoxelIDMap.empty())
     {
         // Diagnostics
         using Triplet    = std::tuple<const recob::Hit*, const recob::Hit*, const recob::Hit*>;
         using TripletMap = std::map<Triplet,std::vector<const recob::SpacePoint*>>;
 
         TripletMap tripletMap;
 
         // So now we loop through the various SpacePoint sources
         for(const auto& spacePointLabel : fSpacePointLabelVec)
         {
             art::Handle< std::vector<recob::SpacePoint>> spacePointHandle;
             event.getByLabel(spacePointLabel, spacePointHandle);
 
             if (!spacePointHandle.isValid()) continue;
 
             // Look up assocations to hits
             art::FindManyP<recob::Hit> spHitAssnVec(spacePointHandle, event, spacePointLabel);
 
             // And now, without further adieu, here we begin the loop that will actually produce some useful output.
             // Loop over all space points and find out their true nature
             for(size_t idx = 0; idx < spacePointHandle->size(); idx++)
             {
                 art::Ptr<recob::SpacePoint> spacePointPtr(spacePointHandle,idx);
 
                 std::vector<art::Ptr<recob::Hit>> associatedHits(spHitAssnVec.at(spacePointPtr.key()));
 
                 if (associatedHits.size() != 3)
                 {
                     mf::LogDebug("SpacePointAnalysisMC") << "I am certain this cannot happen... but here you go, space point with " << associatedHits.size() << " hits" << std::endl;
                     continue;
                 }
 
                 // Retrieve the magic numbers from the space point
                 float spQuality   = spacePointPtr->Chisq();
                 float spCharge    = spacePointPtr->ErrXYZ()[1];
                 float spAsymmetry = spacePointPtr->ErrXYZ()[3];
                 float smallestPH  = std::numeric_limits<float>::max();
                 float largestPH   = 0.;
                 int   numHits     = 0;
                 float averagePH   = 0.;
                 float averagePT   = 0.;
                 float largestDelT = 0.;
 
                 std::vector<int>            numIDEsHitVec;
                 std::vector<int>            simHitDeltaTVec = {0,0,0};
                 std::vector<float>          hitPeakTimeVec  = {-100.,-100.,-100.};
                 std::vector<float>          bigElecDepVec   = {0.,0.,0.};
                 std::vector<unsigned short> bigTDCVec       = {0,0,0};
                 int                         numIDEsSpacePoint(0);
                 int                         numLongHits(0);
                 int                         numIntersections(0);
 
                 std::vector<RecobHitToVoxelIDMap::const_iterator> recobHitToVoxelIterVec;
 
                 std::vector<const recob::Hit*> recobHitVec(3,nullptr);
 
                 // Now we can use our maps to find out if the hits making up the SpacePoint are truly related...
                 for(const auto& hitPtr : associatedHits)
                 {
                     RecobHitToVoxelIDMap::const_iterator hitToVoxelItr = recobHitToVoxelIDMap.find(hitPtr.get());
 
                     float peakAmplitude = hitPtr->PeakAmplitude();
                     float peakTime      = hitPtr->PeakTime();
                     int   plane         = hitPtr->WireID().Plane;
 
                     recobHitVec[plane] = hitPtr.get();
 
                     numHits++;
                     averagePH += peakAmplitude;
                     averagePT += peakTime;
 
                     smallestPH = std::min(peakAmplitude,smallestPH);
                     largestPH  = std::max(peakAmplitude,largestPH);
 
                     if (hitPtr->DegreesOfFreedom() < 2) numLongHits++;
 
                     if (hitToVoxelItr == recobHitToVoxelIDMap.end())
                     {
                         numIDEsHitVec.push_back(0);
                         continue;
                     }
 
                     recobHitToVoxelIterVec.push_back(hitToVoxelItr);
                     numIDEsHitVec.push_back(hitToVoxelItr->second.size());
 
                     hitPeakTimeVec[plane] = clockData.TPCTick2TDC(peakTime);
                 }
 
                 Triplet hitTriplet(recobHitVec[0],recobHitVec[1],recobHitVec[2]);
 
                 tripletMap[hitTriplet].emplace_back(spacePointPtr.get());
 
                 averagePH /= float(numHits);
                 averagePT /= float(numHits);
 
                 for(const auto& hitPtr : associatedHits)
                 {
                     float delT = hitPtr->PeakTime() - averagePT;
 
                     if (std::abs(delT) > std::abs(largestDelT)) largestDelT = delT;
                 }
 
                 // If a SpacePoint is made from "true" MC hits then we will have found the relations to the MC info for all three
                 // hits. If this condition is not satisfied it means one or more hits making the SpacePoint are noise hits
                 if (recobHitToVoxelIterVec.size() == 3)
                 {
                     // Find the intersection of the vectors of IDEs for the first two hits
                     std::vector<sim::LArVoxelID> firstIntersectionVec(recobHitToVoxelIterVec[0]->second.size()+recobHitToVoxelIterVec[1]->second.size());
 
                     std::vector<sim::LArVoxelID>::iterator firstIntersectionItr = std::set_intersection(recobHitToVoxelIterVec[0]->second.begin(),recobHitToVoxelIterVec[0]->second.end(),
                                                                                                         recobHitToVoxelIterVec[1]->second.begin(),recobHitToVoxelIterVec[1]->second.end(),
                                                                                                         firstIntersectionVec.begin());
 
                     firstIntersectionVec.resize(firstIntersectionItr - firstIntersectionVec.begin());
 
                     // No intersection means, of course, the hits did not come from the same MC energy deposit
                     if (!firstIntersectionVec.empty())
                     {
                         // Now find the intersection of the resultant intersection above and the third hit
                         std::vector<sim::LArVoxelID> secondIntersectionVec(firstIntersectionVec.size()+recobHitToVoxelIterVec[2]->second.size());
 
                         std::vector<sim::LArVoxelID>::iterator secondIntersectionItr = std::set_intersection(firstIntersectionVec.begin(),             firstIntersectionVec.end(),
                                                                                                              recobHitToVoxelIterVec[2]->second.begin(),recobHitToVoxelIterVec[2]->second.end(),
                                                                                                              secondIntersectionVec.begin());
 
                         secondIntersectionVec.resize(secondIntersectionItr - secondIntersectionVec.begin());
 
                         numIntersections++;
 
                         // Again, no IDEs in the intersection means it is a ghost space point but, of course, we are hoping
                         // there are common IDEs so we can call it a real SpacePoint
                         if (!secondIntersectionVec.empty())
                         {
                             for(const sim::LArVoxelID& voxelID : secondIntersectionVec)
                             {
                                 VoxelIDToPlaneTDCIDEMap::const_iterator planeToTDCToIDESetMap = voxelToPlaneTDCIDEMap.find(voxelID);
 
                                 if (planeToTDCToIDESetMap->second.size() > 2)
                                 {
                                     numIDEsSpacePoint += 1;
 
                                     for(const auto& planeInfoPair : planeToTDCToIDESetMap->second)
                                     {
                                         unsigned short plane  = planeInfoPair.first;
                                         float          phBig  = 0.;
                                         unsigned short tdcBig = 0;
 
                                         for(const auto& tdcIDEPair : planeInfoPair.second)
                                         {
                                             for(const auto& ide : tdcIDEPair.second)
                                             {
                                                 if (phBig < ide->numElectrons)
                                                 {
                                                     phBig  = ide->numElectrons;
                                                     tdcBig = tdcIDEPair.first;
                                                 }
                                             }
                                         }
 
                                         bigElecDepVec[plane] = phBig;
                                         bigTDCVec[plane]     = tdcBig;
                                     }
                                 }
                                 else std::cout << "   --> Not matching all three planes" << std::endl;
                             }
 
                             numIntersections++;
                         }
                     }
                 }
 
                 // Fill for "all" cases
                 fHitSpacePointObj.fSPQualityVec.push_back(spQuality);
                 fHitSpacePointObj.fSPTotalChargeVec.push_back(spCharge);
                 fHitSpacePointObj.fSPAsymmetryVec.push_back(spAsymmetry);
                 fHitSpacePointObj.fSmallestPHVec.push_back(smallestPH);
                 fHitSpacePointObj.fLargestPHVec.push_back(largestPH);
                 fHitSpacePointObj.fAveragePHVec.push_back(averagePH);
                 fHitSpacePointObj.fLargestDelTVec.push_back(largestDelT);
 
                 fHitSpacePointObj.fNumIDEsHit0Vec.push_back(numIDEsHitVec[0]);
                 fHitSpacePointObj.fNumIDEsHit1Vec.push_back(numIDEsHitVec[1]);
                 fHitSpacePointObj.fNumIDEsHit2Vec.push_back(numIDEsHitVec[2]);
                 fHitSpacePointObj.fNumIDEsSpacePointVec.push_back(numIDEsSpacePoint);
 
                 fHitSpacePointObj.fNumLongHitsVec.emplace_back(numLongHits);
                 fHitSpacePointObj.fNumPlanesSimMatchVec.emplace_back(recobHitToVoxelIterVec.size());
                 fHitSpacePointObj.fNumIntersectSetVec.emplace_back(numIntersections);
                 fHitSpacePointObj.fSimHitDeltaT0Vec.emplace_back(bigTDCVec[0] - hitPeakTimeVec[0]);
                 fHitSpacePointObj.fSimHitDeltaT1Vec.emplace_back(bigTDCVec[1] - hitPeakTimeVec[1]);
                 fHitSpacePointObj.fSimHitDeltaT2Vec.emplace_back(bigTDCVec[2] - hitPeakTimeVec[2]);
 
                 fHitSpacePointObj.fSimDelta10Vec.emplace_back(bigTDCVec[1] - bigTDCVec[0]);
                 fHitSpacePointObj.fSimDelta21Vec.emplace_back(bigTDCVec[2] - bigTDCVec[1]);
                 fHitSpacePointObj.fHitDelta10Vec.emplace_back(hitPeakTimeVec[1] - hitPeakTimeVec[0]);
                 fHitSpacePointObj.fHitDelta21Vec.emplace_back(hitPeakTimeVec[2] - hitPeakTimeVec[1]);
                 fHitSpacePointObj.fBigElecDep0Vec.emplace_back(bigElecDepVec[0]);
                 fHitSpacePointObj.fBigElecDep1Vec.emplace_back(bigElecDepVec[1]);
                 fHitSpacePointObj.fBigElecDep2Vec.emplace_back(bigElecDepVec[2]);
             }
         }
 
         // Can we check to see if we have duplicates?
         std::vector<int> numSpacePointVec = {0,0,0,0,0};
         for(const auto& tripletPair : tripletMap)
         {
             int numSpacePoints = std::min(numSpacePointVec.size()-1,tripletPair.second.size());
             numSpacePointVec[numSpacePoints]++;
         }
         std::cout << "====>> Found " << tripletMap.size() << " SpacePoints, numbers: ";
         for(const auto& count : numSpacePointVec) std::cout << count << " ";
         std::cout << std::endl;
     }
 
     return;
 }

void SpacePointAnalysisMC::SpacePointAnalysisMC::configure ( fhicl::ParameterSet const & pset )

overridevirtual

Reconfigure method.

Arguments:

pset - Fcl parameter set.

Implements IHitEfficiencyHistogramTool.

Definition at line 464 of file SpacePointAnalysisMC_tool.cc.

 {
     fRecobHitLabelVec         = pset.get< std::vector<art::InputTag>>("HitLabelVec",         std::vector<art::InputTag>() = {"cluster3d"});
     fSpacePointLabelVec       = pset.get< std::vector<art::InputTag>>("SpacePointLabelVec",  std::vector<art::InputTag>() = {"cluster3d"});
     fMCParticleProducerLabel  = pset.get< art::InputTag             >("MCParticleLabel",     "largeant");
     fSimChannelProducerLabel  = pset.get< art::InputTag             >("SimChannelLabel",     "largeant");
     fSimEnergyProducerLabel   = pset.get< art::InputTag             >("SimEnergyLabel",      "largeant");
     fOffsetVec                = pset.get<std::vector<int>           >("OffsetVec",           std::vector<int>()={0,0,0});
     fSimChannelMinEnergy      = pset.get<float                      >("SimChannelMinEnergy", std::numeric_limits<float>::epsilon());
     fSimEnergyMinEnergy       = pset.get<float                      >("SimEnergyMinEnergy",  std::numeric_limits<float>::epsilon());
 
     return;
 }

void SpacePointAnalysisMC::SpacePointAnalysisMC::endJob ( int numEvents )

overridevirtual

Interface for method to executve at the end of run processing.

Parameters

int	number of events to use for normalization

Implements IHitEfficiencyHistogramTool.

Definition at line 1186 of file SpacePointAnalysisMC_tool.cc.

 {
     return;
 }

void SpacePointAnalysisMC::SpacePointAnalysisMC::fillHistograms ( const art::Event & event ) const

overridevirtual

Interface for filling histograms.

Implements IHitEfficiencyHistogramTool.

Definition at line 533 of file SpacePointAnalysisMC_tool.cc.

 {
     // Ok... this is starting to grow too much and get out of control... we will need to break it up directly...
 
     // Always clear the tuple
     clear();
 
     art::Handle<std::vector<sim::SimChannel>> simChannelHandle;
     event.getByLabel(fSimChannelProducerLabel, simChannelHandle);
 
     if (!simChannelHandle.isValid() || simChannelHandle->empty() ) return;
 
     art::Handle<std::vector<sim::SimEnergyDeposit>> simEnergyHandle;
     event.getByLabel(fSimEnergyProducerLabel, simEnergyHandle);
 
 //    if (!simEnergyHandle.isValid() || simEnergyHandle->empty()) return;
 
     art::Handle<std::vector<simb::MCParticle>> mcParticleHandle;
     event.getByLabel(fMCParticleProducerLabel, mcParticleHandle);
 
     // If there is no sim channel informaton then exit
     if (!mcParticleHandle.isValid()) return;
 
     mf::LogDebug("SpacePointAnalysisMC") << "+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++" << std::endl;
 
     // First task is to build a map between ides and voxel ids (that we calcualate based on position)
     // and also get the reverse since it will be useful in the end.
     // At the same time should also build a mapping of ides per channel so we can do quick hit lookup
     IDEToVoxelIDMap         ideToVoxelIDMap;
     VoxelIDToIDESetMap      voxelIDToIDEMap;
     ChanToTDCToIDEMap       chanToTDCToIDEMap;
     VoxelIDSet              simChannelVoxelIDSet;
     VoxelIDToPlaneTDCIDEMap voxelIDToPlaneTDCIDEMap;
 
     TrackIDChanToTDCIDEMap trackIDChanToTDCIDEMap;
 
     // Fill the above maps/structures
     for(const auto& simChannel : *simChannelHandle)
     {
         raw::ChannelID_t channel = simChannel.Channel();
 
         geo::WireID wireID = fGeometry->ChannelToWire(channel).front();
 
         for(const auto& tdcide : simChannel.TDCIDEMap())
         {
             for(const auto& ide : tdcide.second) //chanToTDCToIDEMap[simChannel.Channel()][tdcide.first] = ide;
             {
                 if (ide.energy < fSimChannelMinEnergy) continue;
 
                 sim::LArVoxelID voxelID(ide.x,ide.y,ide.z,0.);
 
                 ideToVoxelIDMap[&ide]    = voxelID;
                 voxelIDToIDEMap[voxelID].insert(&ide);
                 chanToTDCToIDEMap[simChannel.Channel()][tdcide.first].insert(&ide);
                 simChannelVoxelIDSet.insert(voxelID);
 
                 trackIDChanToTDCIDEMap[ide.trackID][simChannel.Channel()].emplace_back(tdcide.first,&ide);
 
                 voxelIDToPlaneTDCIDEMap[voxelID][wireID.Plane][tdcide.first].insert(&ide);
 
                 if (ide.energy < std::numeric_limits<float>::epsilon()) mf::LogDebug("SpacePointAnalysisMC") << ">> epsilon simchan deposited energy: " << ide.energy << std::endl;
             }
         }
     }
 
     // More data structures, here we want to keep track of the start/peak/end of the charge deposit along a wire for a given track
     TrackToChanChargeMap trackToChanChargeMap;
 
     // Go through the list of track to ides and get the total deposited energy per track
     float bestTotDepEne(0.);
     int   bestTrackID(0);
 
     makeTrackToChanChargeMap(trackIDChanToTDCIDEMap, trackToChanChargeMap, bestTotDepEne, bestTrackID);
 
     // Ok, for my next trick I want to build a mapping between hits and voxel IDs. Note that any given hit can be associated to more than one voxel...
     // We do this on the entire hit collection, ultimately we will want to consider SpacePoint efficiency (this could be done in the loop over SpacePoints
     // using the associated hits and would save time/memory)
     using VoxelIDSet           = std::set<sim::LArVoxelID>;
     using RecobHitToVoxelIDMap = std::unordered_map<const recob::Hit*, VoxelIDSet>;
 
     RecobHitToVoxelIDMap recobHitToVoxelIDMap;
 
     ChanToTDCIDEMap& chanToTDCIDEMap = trackIDChanToTDCIDEMap[bestTrackID];
 
     // Recover the "best" track info to start
     TrackToChanChargeMap::const_iterator chanToChargeMapItr = trackToChanChargeMap.find(bestTrackID);
 
     // Process the hit/simulation
     compareHitsToSim(event, chanToTDCToIDEMap, chanToChargeMapItr->second, chanToTDCIDEMap, ideToVoxelIDMap, recobHitToVoxelIDMap);
 
     // Now do the space points
     auto const clockData = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(event);
     compareSpacePointsToSim(event, clockData, recobHitToVoxelIDMap, voxelIDToPlaneTDCIDEMap);
 
     // Make sure the output tuples are filled
     fHitSpacePointObj.fill();
 
     for(auto& hitObj : fHitSimObjVec) hitObj.fill();
 
     return;
 }

void SpacePointAnalysisMC::SpacePointAnalysisMC::initializeHists	(	art::ServiceHandle< art::TFileService > &	tfs,
		const std::string &	dirName
	)

overridevirtual

Interface for initializing the histograms to be filled.

Begin job method.

Parameters

TFileService	handle to the TFile service
string	subdirectory to store the hists in

Implements IHitEfficiencyHistogramTool.

Definition at line 480 of file SpacePointAnalysisMC_tool.cc.

 {
     // Make a directory for these histograms
 //    art::TFileDirectory dir = tfs->mkdir(dirName.c_str());
 
     return;
 }

void SpacePointAnalysisMC::SpacePointAnalysisMC::initializeTuple ( TTree * tree )

overridevirtual

Interface for initializing the tuple variables.

Parameters

TTree pointer to a TTree object to which to add variables

Implements IHitEfficiencyHistogramTool.

Definition at line 488 of file SpacePointAnalysisMC_tool.cc.

 {
     // Access ART's TFileService, which will handle creating and writing
     // histograms and n-tuples for us.
     art::ServiceHandle<art::TFileService> tfs;
 
     fTree = tree;
 
     fTree->Branch("CryostataVec",       "std::vector<int>",   &fCryoVec);
     fTree->Branch("TPCVec",             "std::vector<int>",   &fTPCVec);
     fTree->Branch("PlaneVec",           "std::vector<int>",   &fPlaneVec);
 
     // Set up specific branch for space points
     TTree* locTree = tfs->makeAndRegister<TTree>("SpacePoint_t","SpacePoint Tuple");
 
     fHitSpacePointObj.setBranches(locTree);
 
     fHitSimObjVec.resize(fGeometry->Nplanes());
 
     for(size_t plane = 0; plane < fGeometry->Nplanes(); plane++)
     {
         // Set up specific branch for space points
         locTree = tfs->makeAndRegister<TTree>("MatchedHits_P"+std::to_string(plane),"Matched Hits Tuple plane "+std::to_string(plane));
 
         fHitSimObjVec[plane].setBranches(locTree);
     }
 
     clear();
 
     return;
 }

void SpacePointAnalysisMC::SpacePointAnalysisMC::makeTrackToChanChargeMap	(	const TrackIDChanToTDCIDEMap &	trackIDChanToTDCIDEMap,
		TrackToChanChargeMap &	trackToChanChargeMap,
		float &	bestTotDepEne,
		int &	bestTrackID
	)		const

private

Definition at line 635 of file SpacePointAnalysisMC_tool.cc.

 {
     // Pretty straightforward looping here...
     for(const auto& trackIDEPair : trackIDChanToTDCIDEMap)
     {
         ChanToChargeMap& chanToChargeMap = trackToChanChargeMap[trackIDEPair.first];
 
         float trackTotDepE(0.);
 
         for(const auto& chanTDCIDEPair : trackIDEPair.second)
         {
             ChargeDepositVec& chargeDepositVec = chanToChargeMap[chanTDCIDEPair.first];
 
             // Keep track of first,peak,last/ene
             TDCIDEPair firstPair = chanTDCIDEPair.second.front();
             TDCIDEPair peakPair  = firstPair;
             TDCIDEPair lastPair  = chanTDCIDEPair.second.back();
 
             // Keep watch for gaps
             TDCIDEPair prevPair  = firstPair;
 
             // Keep track of deposited energy on a snippet
             float snippetDepEne(0.);
             float snippetNumElectrons(0.);
 
             for(const auto& tdcIDEPair : chanTDCIDEPair.second)
             {
                 float depEne = tdcIDEPair.second->energy;
 
                 trackTotDepE += depEne;
 
                 // Watch for a gap...
                 if (tdcIDEPair.first - prevPair.first > 1)
                 {
                     chargeDepositVec.emplace_back(firstPair,peakPair,prevPair,snippetDepEne,snippetNumElectrons);
 
                     firstPair           = tdcIDEPair;
                     peakPair            = firstPair;
                     snippetDepEne       = 0.;
                     snippetNumElectrons = 0.;
                 }
 
                 snippetDepEne       += depEne;
                 snippetNumElectrons += tdcIDEPair.second->numElectrons;
 
                 if (depEne > peakPair.second->energy) peakPair = tdcIDEPair;
 
                 prevPair = tdcIDEPair;
             }
 
             chargeDepositVec.emplace_back(firstPair,peakPair,lastPair,snippetDepEne,snippetNumElectrons);
         }
 
         if (trackTotDepE > bestTotDepEne)
         {
             bestTrackID   = trackIDEPair.first;
             bestTotDepEne = trackTotDepE;
         }
     }
 
     return;
 }

void SpacePointAnalysisMC::SpacePointAnalysisMC::matchHitSim	(	const detinfo::DetectorClocksData &	clockData,
		const HitPointerVec &	hitPointerVec,
		const ChanToTDCToIDEMap &	chanToTDCToIDEMap,
		const ChargeDepositVec &	chargeDepositVec,
		const ChanToTDCIDEMap &	chanToTDCIDEMap,
		const IDEToVoxelIDMap &	ideToVoxelIDMap,
		RecobHitToVoxelIDMap &	recobHitToVoxelIDMap
	)		const

private

Definition at line 794 of file SpacePointAnalysisMC_tool.cc.

 {
     // Data structure to allow ordering of multiple hits in a snippet
     using HitPeakTimeChargeTuple = std::tuple<int,const recob::Hit*,ChargeDepositVec::const_iterator>;
     using HitPeakTimeChargeVec   = std::vector<HitPeakTimeChargeTuple>;
 
     HitPeakTimeChargeVec hitPeakTimeChargeVec;
     float                maxPulseHeight(1.);   // Don't let this be zero
 
     // If here then we are on to the next hit, so we need to process our current list
     for(const auto& hit : hitPointerVec)
     {
         // Recover hit time range (in ticks), cast a wide net here
         int peakTick  = std::round(hit->PeakTime());
         int startTick = std::max(   0,int(std::floor(hit->PeakTime() - 3. * hit->RMS())));
         int endTick   = std::min(4096,int(std::ceil(hit->PeakTime() + 3. * hit->RMS())));
 
         int startTDC = clockData.TPCTick2TDC(startTick - fOffsetVec[hit->WireID().Plane]);
         int peakTDC  = clockData.TPCTick2TDC(peakTick  - fOffsetVec[hit->WireID().Plane]);
         int endTDC   = clockData.TPCTick2TDC(endTick   - fOffsetVec[hit->WireID().Plane]);
 
         maxPulseHeight = std::max(maxPulseHeight,hit->PeakAmplitude());
 
         // If we have a match then this iterator gets set to the matching values
         ChargeDepositVec::const_iterator chargeMatchItr = chargeDepositVec.end();
 
         int   bestPeakDiff = std::numeric_limits<int>::max();
 
         // Match the hit (if there is one)
         for(ChargeDepositVec::const_iterator chargeInfoItr = chargeDepositVec.begin(); chargeInfoItr != chargeDepositVec.end(); chargeInfoItr++)
         {
             // Require some amount of overlap between the hit and the sim info
             if (endTDC > std::get<0>(*chargeInfoItr).first && startTDC < std::get<2>(*chargeInfoItr).first)
             {
                 const TDCIDEPair& peakTDCIDE = std::get<1>(*chargeInfoItr);
 
                 int peakDiff = peakTDC - int(peakTDCIDE.first);
 
                 if (std::abs(peakDiff) < std::abs(bestPeakDiff))
                 {
                     bestPeakDiff   = peakDiff;
                     chargeMatchItr = chargeInfoItr;
                 }
             }
         }
 
         // If no match then skip
         if (chargeMatchItr == chargeDepositVec.end()) continue;
 
         hitPeakTimeChargeVec.emplace_back(std::make_tuple(bestPeakDiff,hit,chargeMatchItr));
     }
 
     if (!hitPeakTimeChargeVec.empty())
     {
         // Ok, now we sort this vector by smallest peak time
         std::sort(hitPeakTimeChargeVec.begin(),hitPeakTimeChargeVec.end(),[](const auto& left,const auto& right){return std::abs(std::get<0>(left)) < std::abs(std::get<0>(right));});
 
         // Keep track of hit ordering on this snippet
         int hitOrder(0);
 
         HitSimulationTupleObj& hitObj = fHitSimObjVec[std::get<1>(hitPeakTimeChargeVec.front())->WireID().Plane];
 
         // Now loop through
         for(const auto& hitPeakCharge : hitPeakTimeChargeVec)
         {
             const recob::Hit*    hit           = std::get<1>(hitPeakCharge);
             const ChargeDeposit& chargeDeposit = *std::get<2>(hitPeakCharge);
 
             // Recover hit time range (in ticks), cast a wide net here
             int   peakTick  = std::round(hit->PeakTime());
             int   startTick = std::max(   0,int(std::floor(hit->PeakTime() - 3. * hit->RMS())));
             int   endTick   = std::min(4096,int(std::ceil(hit->PeakTime() + 3. * hit->RMS())));
 
             int   startTDC = clockData.TPCTick2TDC(startTick - fOffsetVec[hit->WireID().Plane]);
             int   peakTDC  = clockData.TPCTick2TDC(peakTick  - fOffsetVec[hit->WireID().Plane]);
             int   endTDC   = clockData.TPCTick2TDC(endTick   - fOffsetVec[hit->WireID().Plane]);
 
             int   firstSimTick(std::get<0>(chargeDeposit).first);
             int   lastSimTick(std::get<2>(chargeDeposit).first);
             int   maxDepTick(std::get<1>(chargeDeposit).first);
             float maxDepEneTick(std::get<1>(chargeDeposit).second->energy);
             float bestNumElectrons(std::get<4>(chargeDeposit));
             float bestDepEne(std::get<3>(chargeDeposit));
             float totDepEne(0.);
             float totNumElectrons(0.);
             int   bestTicks(lastSimTick - firstSimTick + 1);
 
             // We want to get the total energy deposit from all particles in the ticks for this hit
             const TDCToIDEMap& tdcToIDEMap = chanToTDCToIDEMap.find(hit->Channel())->second;
             for(const auto& tdcToIDEPair : tdcToIDEMap)
             {
                 for(const auto& ide : tdcToIDEPair.second)
                 {
                     totDepEne       += ide->energy;
                     totNumElectrons += ide->numElectrons;
                 }
             }
 
             // One final time through to find sim ticks that "matter"
             // We define this as the collection of IDE's that make up to 90% of the total deposit
             ChanToTDCIDEMap::const_iterator tickToTDCIDEVecItr = chanToTDCIDEMap.find(hit->Channel());
 
             if (tickToTDCIDEVecItr == chanToTDCIDEMap.end()) continue;
 
             const TickTDCIDEVec& tickToTDCIDEVec = tickToTDCIDEVecItr->second;
             TickTDCIDEVec        tickIDEVec;
 
             for(const auto& tickInfo : tickToTDCIDEVec)
             {
                 //if (tickInfo.first >= firstSimTick && tickInfo.first <= lastSimTick) tickIDEVec.emplace_back(tickInfo);
                 if (tickInfo.first >= startTDC && tickInfo.first <= endTDC) tickIDEVec.emplace_back(tickInfo);
             }
 
             std::sort(tickIDEVec.begin(),tickIDEVec.end(),[](const auto& left,const auto& right){return left.second->energy > right.second->energy;});
 
             // Grab the voxelID set for this tick
             VoxelIDSet voxelIDSet;
 
             int   bestTicksGood(0);
             float sumEne(0.);
 
             for(const auto& tickInfo : tickIDEVec)
             {
                 // At the same time we can keep track of the voxels associated to the best track
                 IDEToVoxelIDMap::const_iterator ideToVoxelIDMapItr = ideToVoxelIDMap.find(tickInfo.second);
 
                 if (ideToVoxelIDMapItr == ideToVoxelIDMap.end()) continue;
 
                 const sim::LArVoxelID& voxelID = ideToVoxelIDMapItr->second;
 
                 sumEne += tickInfo.second->energy;
                 bestTicksGood++;
 
                 voxelIDSet.insert(voxelID);
 
                 if (sumEne > 0.9 * bestDepEne) break;
             }
 
             // Finally, grab the voxels from the track leaving the most energy
             recobHitToVoxelIDMap[hit] = voxelIDSet;
 
             hitObj.fillSimInfo(bestTicks, bestTicksGood, totDepEne, totNumElectrons, bestDepEne, bestNumElectrons, maxDepEneTick);
             hitObj.fillMixedInfo(endTick-startTick+1, maxDepTick-startTDC, peakTDC-maxDepTick);
             hitObj.fillHitInfo(hit,hitOrder++);
         }
 
         // Resort in pulse height order (largest to smallest)
         std::sort(hitPeakTimeChargeVec.begin(),hitPeakTimeChargeVec.end(),[](const auto& left,const auto& right){return std::get<1>(left)->PeakAmplitude() > std::get<1>(right)->PeakAmplitude();});
 
         // Now loop through
         for(const auto& hitPeakCharge : hitPeakTimeChargeVec)
         {
             hitObj.fPHOrderHitVec.emplace_back(std::get<1>(hitPeakCharge)->LocalIndex());
             hitObj.fPHFractionVec.emplace_back(std::get<1>(hitPeakCharge)->PeakAmplitude() / maxPulseHeight);
         }
     }
 
     return;
 }