Inheritance diagram for SpacePointAnalysis::SpacePointAnalysis:

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

	~SpacePointAnalysis ()
	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	HitPointerVec = std::vector< const recob::Hit * >

using	PlaneToT0OffsetMap = std::map< geo::PlaneID, float >

using	HitTuplebjVec = std::vector< HitTupleObj >

Private Member Functions
void	clear () const

void	processSpacePoints (const art::Event &, const detinfo::DetectorClocksData &) const

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

art::InputTag	fBadChannelProducerLabel

bool	fUseT0Offsets

TTree *	fTree

std::vector< int >	fTPCVec

std::vector< int >	fCryoVec

std::vector< int >	fPlaneVec

PlaneToT0OffsetMap	fPlaneToT0OffsetMap

HitTuplebjVec	fHitTupleObjVec

HitSpacePointObj	fHitSpacePointObj

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

Detailed Description

Definition at line 264 of file SpacePointAnalysis_tool.cc.

Member Typedef Documentation

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

private

Definition at line 317 of file SpacePointAnalysis_tool.cc.

using SpacePointAnalysis::SpacePointAnalysis::HitTuplebjVec = std::vector<HitTupleObj>

private

Definition at line 335 of file SpacePointAnalysis_tool.cc.

using SpacePointAnalysis::SpacePointAnalysis::PlaneToT0OffsetMap = std::map<geo::PlaneID,float>

private

Definition at line 320 of file SpacePointAnalysis_tool.cc.

Constructor & Destructor Documentation

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

explicit

Constructor.

Parameters

pset	Constructor.

Arguments:

pset - Fcl parameters.

Definition at line 351 of file SpacePointAnalysis_tool.cc.

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

SpacePointAnalysis::SpacePointAnalysis::~SpacePointAnalysis ( )

Destructor.

Definition at line 363 of file SpacePointAnalysis_tool.cc.

364 {}

Member Function Documentation

void SpacePointAnalysis::SpacePointAnalysis::clear ( ) const

private

Definition at line 423 of file SpacePointAnalysis_tool.cc.

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

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

overridevirtual

Reconfigure method.

Arguments:

pset - Fcl parameter set.

Implements IHitEfficiencyHistogramTool.

Definition at line 373 of file SpacePointAnalysis_tool.cc.

 {
     fSpacePointLabelVec  = pset.get< std::vector<art::InputTag>>("SpacePointLabelVec",  {"cluster3d"});
     fUseT0Offsets        = pset.get< bool                      >("UseT0Offsets",        false);
 
     return;
 }

void SpacePointAnalysis::SpacePointAnalysis::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 690 of file SpacePointAnalysis_tool.cc.

 {
     return;
 }

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

overridevirtual

Interface for filling histograms.

Implements IHitEfficiencyHistogramTool.

Definition at line 436 of file SpacePointAnalysis_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();
 
     mf::LogDebug("SpacePointAnalysis") << "+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++" << std::endl;
 
     // Now do the space points
     auto const clockData = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(event);
     processSpacePoints(event, clockData);
 
     // Make sure the output tuples are filled
     fHitSpacePointObj.fill();
 
     for(auto& hitObj : fHitTupleObjVec) hitObj.fill();
 
     return;
 }

void SpacePointAnalysis::SpacePointAnalysis::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 383 of file SpacePointAnalysis_tool.cc.

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

void SpacePointAnalysis::SpacePointAnalysis::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 391 of file SpacePointAnalysis_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);
 
     fHitTupleObjVec.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));
 
         fHitTupleObjVec[plane].setBranches(locTree);
     }
 
     clear();
 
     return;
 }

void SpacePointAnalysis::SpacePointAnalysis::processSpacePoints	(	const art::Event &	event,
		const detinfo::DetectorClocksData &	clockData
	)		const

private

Definition at line 457 of file SpacePointAnalysis_tool.cc.

 {
     // One time initialization done?
 
     // Do the one time initialization of the tick offsets. 
     if (fPlaneToT0OffsetMap.empty())
     {
         // Need the detector properties which needs the clocks 
         auto const det_prop = art::ServiceHandle<detinfo::DetectorPropertiesService const>()->DataFor(event, clockData);
 
         for(size_t cryoIdx = 0; cryoIdx < fGeometry->Ncryostats(); cryoIdx++)
         {
             for(size_t tpcIdx = 0; tpcIdx < fGeometry->NTPC(); tpcIdx++)
             {
                 for(size_t planeIdx = 0; planeIdx < fGeometry->Nplanes(); planeIdx++)
                 {
                     geo::PlaneID planeID(cryoIdx,tpcIdx,planeIdx);
 
 //                    if (fUseT0Offsets) fPlaneToT0OffsetMap[planeID] = int(det_prop.GetXTicksOffset(planeID) - det_prop.GetXTicksOffset(geo::PlaneID(cryoIdx,tpcIdx,0)));
                     if (fUseT0Offsets) fPlaneToT0OffsetMap[planeID] = det_prop.GetXTicksOffset(planeID) - det_prop.GetXTicksOffset(geo::PlaneID(cryoIdx,tpcIdx,0));
                     else               fPlaneToT0OffsetMap[planeID] = 0.;
 
                     std::cout << "--PlaneID: " << planeID << ", has T0 offset: " << fPlaneToT0OffsetMap.find(planeID)->second << std::endl;
                 }
             }
         }   
     }
 
     // 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;
 
     // detector properties
     auto const detProp = art::ServiceHandle<detinfo::DetectorPropertiesService const>()->DataForJob();
 
     // So now we loop through the various SpacePoint sources
     for(const auto& collectionLabel : fSpacePointLabelVec)
     {
         art::Handle<std::vector<recob::PFParticle>> pfParticleHandle;
         event.getByLabel(collectionLabel, pfParticleHandle);
 
         if (!pfParticleHandle.isValid()) continue;
 
         art::Handle< std::vector<recob::SpacePoint>> spacePointHandle;
         event.getByLabel(collectionLabel, spacePointHandle);
 
         if (!spacePointHandle.isValid()) continue;
 
         // Recover the collection of associations between tracks and hits
         art::FindManyP<recob::SpacePoint> pfPartSpacePointAssns(pfParticleHandle, event, collectionLabel);
 
         // Look up assocations between pfparticles and space points
         art::FindManyP<recob::Hit> spHitAssnVec(spacePointHandle, event, collectionLabel);
 
         // Use this to build a map between PFParticles and the number of associated space points 
         using PFParticleToNumSPMap = std::map<const recob::PFParticle*,int>;
 
         PFParticleToNumSPMap pfParticleToNumSPMap;
 
         for(size_t idx = 0; idx < pfParticleHandle->size(); idx++)
         {
             art::Ptr<recob::PFParticle> pfParticle(pfParticleHandle,idx);
 
             std::vector<art::Ptr<recob::SpacePoint>> spacePointVec(pfPartSpacePointAssns.at(pfParticle.key()));
 
             pfParticleToNumSPMap[pfParticle.get()] = spacePointVec.size();
 
 //            if (spacePointVec.size() > 50000) std::cout << "SpacePointAnalysis finds PFParticle with " << spacePointVec.size() << " associated space points, run/event: " << event.id() << std::endl;
         }
 
         // Ok now we want the reverse look up
         art::FindManyP<recob::PFParticle> spacePointPFPartAssns(spacePointHandle, event, collectionLabel);
 
         std::unordered_map<const recob::Hit*,int> uniqueHitMap;
 
         // 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++)
         {
             // Recover space point
             art::Ptr<recob::SpacePoint> spacePointPtr(spacePointHandle,idx);
 
             std::vector<art::Ptr<recob::Hit>> associatedHits(spHitAssnVec.at(spacePointPtr.key()));
 
             if (associatedHits.size() < 2)
             {
                 mf::LogDebug("SpacePointAnalysis") << "I am certain this cannot happen... but here you go, space point with " << associatedHits.size() << " hits" << std::endl;
                 continue;
             }
 
             // Recover the PFParticle associated to this space point and get the number of associated hits
             int nSpacePointsInPFParticle(0);
 
             std::vector<art::Ptr<recob::PFParticle>> pfParticleVec(spacePointPFPartAssns.at(spacePointPtr.key()));
 
             if (pfParticleVec.size() == 1) nSpacePointsInPFParticle = pfParticleToNumSPMap.at(pfParticleVec[0].get());
 
             mf::LogDebug("SpacePointAnalysis") << "==> pfPartVec size: " << pfParticleVec.size() << ", # space points: " << nSpacePointsInPFParticle << std::endl;
 
             // Retrieve the magic numbers from the space point
             float              spQuality       = spacePointPtr->Chisq();
             float              spCharge        = spacePointPtr->ErrXYZ()[1];
             float              spAsymmetry     = spacePointPtr->ErrXYZ()[3];
             const Double32_t*  spPosition      = spacePointPtr->XYZ();
             float              smallestPH      = std::numeric_limits<float>::max();
             float              largestPH       = 0.;
             int                numHits         = 0;
             float              averagePH       = 0.;
             float              averagePT       = 0.;
             float              largestDelT     = 0.;
             float              smallestDelT    = 100000.;
             std::vector<float> hitPeakTimeVec  = {-10000.,-20000.,-30000.};
             std::vector<float> hitPeakRMSVec   = {1000.,1000.,1000.};
             int                hitMultProduct  = 1;
             int                numLongHits(0);
             int                numIntersections(0);
             int                cryostat(-1);
             int                tpc(-1);
 
             std::vector<const recob::Hit*> recobHitVec(3,nullptr);
 
             std::vector<float> peakAmpVec;
 
             // Now we can use our maps to find out if the hits making up the SpacePoint are truly related...
             for(const auto& hitPtr : associatedHits)
             {
                 float peakAmplitude = hitPtr->PeakAmplitude();
                 float peakTime      = hitPtr->PeakTime();
                 float rms           = hitPtr->RMS();
                 int   plane         = hitPtr->WireID().Plane;
 
                 peakAmpVec.emplace_back(peakAmplitude);
 
                 // Add to the set
                 uniqueHitMap[hitPtr.get()] = nSpacePointsInPFParticle;
 
                 recobHitVec[plane] = hitPtr.get();
                 numHits++;
                 averagePH  += peakAmplitude;
                 smallestPH  = std::min(peakAmplitude,smallestPH);
                 largestPH   = std::max(peakAmplitude,largestPH);
 
                 hitMultProduct *= hitPtr->Multiplicity();
 
                 if (hitPtr->DegreesOfFreedom() < 2) numLongHits++;
 
                 hitPeakTimeVec[plane] = peakTime - fPlaneToT0OffsetMap.find(hitPtr->WireID().planeID())->second;
                 hitPeakRMSVec[plane]  = rms;
                 averagePT            += hitPeakTimeVec[plane];
 
                 cryostat = hitPtr->WireID().Cryostat;
                 tpc      = hitPtr->WireID().TPC;
             }
 
             Triplet hitTriplet(recobHitVec[0],recobHitVec[1],recobHitVec[2]);
 
             tripletMap[hitTriplet].emplace_back(spacePointPtr.get());
 
             averagePH /= float(numHits);
             averagePT /= float(numHits);
 
             for(size_t planeIdx = 0; planeIdx < 3; planeIdx++)
             {
                 // Skip if hit missing
                 if (hitPeakTimeVec[planeIdx] < 0) continue;
 
                 float delT = hitPeakTimeVec[planeIdx] - averagePT;
                 if (std::abs(delT) > std::abs(largestDelT))  largestDelT  = delT;
                 if (std::abs(delT) < std::abs(smallestDelT)) smallestDelT = delT;
             }
 
             // Fill for "all" cases
             fHitSpacePointObj.fSPCryostatVec.emplace_back(cryostat);
             fHitSpacePointObj.fSPTPCVec.emplace_back(tpc);
             fHitSpacePointObj.fSPQualityVec.emplace_back(spQuality);
             fHitSpacePointObj.fSPTotalChargeVec.emplace_back(spCharge);
             fHitSpacePointObj.fSPAsymmetryVec.emplace_back(spAsymmetry);
             fHitSpacePointObj.fSmallestPHVec.emplace_back(smallestPH);
             fHitSpacePointObj.fLargestPHVec.emplace_back(largestPH);
             fHitSpacePointObj.fAveragePHVec.emplace_back(averagePH);
             fHitSpacePointObj.fLargestDelTVec.emplace_back(largestDelT);
             fHitSpacePointObj.fSmallestDelTVec.emplace_back(smallestDelT);
             fHitSpacePointObj.fNum2DHitsVec.emplace_back(numHits);
             fHitSpacePointObj.fNumLongHitsVec.emplace_back(numLongHits);
             fHitSpacePointObj.fNumIntersectSetVec.emplace_back(numIntersections);
             fHitSpacePointObj.fClusterNSPVec.emplace_back(nSpacePointsInPFParticle);
             fHitSpacePointObj.fHitDelta10Vec.emplace_back(hitPeakTimeVec[1] - hitPeakTimeVec[0]);
             fHitSpacePointObj.fHitSigma10Vec.emplace_back(std::sqrt(std::pow(hitPeakRMSVec[1],2) + std::pow(hitPeakRMSVec[0],2)));
             fHitSpacePointObj.fHitDelta21Vec.emplace_back(hitPeakTimeVec[2] - hitPeakTimeVec[1]);
             fHitSpacePointObj.fHitSigma21Vec.emplace_back(std::sqrt(std::pow(hitPeakRMSVec[2],2) + std::pow(hitPeakRMSVec[1],2)));
             fHitSpacePointObj.fHitDelta20Vec.emplace_back(hitPeakTimeVec[2] - hitPeakTimeVec[0]);
             fHitSpacePointObj.fHitSigma20Vec.emplace_back(std::sqrt(std::pow(hitPeakRMSVec[2],2) + std::pow(hitPeakRMSVec[0],2)));
             fHitSpacePointObj.fHitMultProductVec.emplace_back(hitMultProduct);
 
             fHitSpacePointObj.fSP_x.emplace_back(spPosition[0]);
             fHitSpacePointObj.fSP_y.emplace_back(spPosition[1]);
             fHitSpacePointObj.fSP_z.emplace_back(spPosition[2]);
         }
 
         // Now include hit information for unique hits
         for(const auto& hitItr : uniqueHitMap)
         {
             // Recover hit time range (in ticks), cast a wide net here
             const recob::Hit* hit = hitItr.first;
 
             float peakTime  = hit->PeakTime();
             float rms       = hit->PeakTime();
             int   startTick = std::max(   0,int(std::floor(peakTime - 3. * rms)));
             int   endTick   = std::min(4096,int(std::ceil( peakTime + 3. * rms)));
 
             fHitTupleObjVec[hit->WireID().Plane].fillHitInfo(hit,endTick-startTick+1,hitItr.second);
 
         }
     }
     // 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;
 }