3 #include "fhiclcpp/ParameterSet.h"
4 #include "art/Utilities/ToolMacros.h"
5 #include "art/Framework/Services/Registry/ServiceHandle.h"
6 #include "art_root_io/TFileService.h"
7 #include "art/Framework/Core/ModuleMacros.h"
8 #include "art_root_io/TFileDirectory.h"
9 #include "messagefacility/MessageLogger/MessageLogger.h"
21 #include "nusimdata/SimulationBase/MCParticle.h"
24 #include <Eigen/Dense>
29 #include "TProfile2D.h"
35 namespace TrackHitEfficiencyAnalysis
55 using HitPtrVec = std::vector<art::Ptr<recob::Hit>>;
75 void configure(fhicl::ParameterSet
const & pset)
override;
83 void initializeHists(art::ServiceHandle<art::TFileService>&,
const std::string&)
override;
97 void endJob(
int numEvents)
override;
157 mutable TTree*
fTree;
159 mutable std::vector<int>
fTPCVec;
202 fGeometry = lar::providerFrom<geo::Geometry>();
207 mf::LogInfo(
"TrackHitEfficiencyAnalysis") <<
"TrackHitEfficiencyAnalysis configured\n";
212 TrackHitEfficiencyAnalysis::~TrackHitEfficiencyAnalysis()
222 void TrackHitEfficiencyAnalysis::configure(fhicl::ParameterSet
const & pset)
224 fRawDigitProducerLabelVec = pset.get< std::vector<art::InputTag>>(
"RawDigitLabelVec", std::vector<art::InputTag>() = {
"rawdigitfilter"});
225 fWireProducerLabelVec = pset.get< std::vector<art::InputTag>>(
"WireModuleLabelVec", std::vector<art::InputTag>() = {
"decon1droi"});
226 fHitProducerLabelVec = pset.get< std::vector<art::InputTag>>(
"HitModuleLabelVec", std::vector<art::InputTag>() = {
"gauss"});
231 fLocalDirName = pset.get<std::string >(
"LocalDirName", std::string(
"wow"));
232 fOffsetVec = pset.get<std::vector<int> >(
"OffsetVec", std::vector<int>()={0,0,0});
233 fSigmaVec = pset.get<std::vector<float> >(
"SigmaVec", std::vector<float>()={1.,1.,1.});
239 void TrackHitEfficiencyAnalysis::initializeHists(art::ServiceHandle<art::TFileService>&
tfs,
const std::string& dirName)
242 art::TFileDirectory
dir = tfs->mkdir(dirName.c_str());
297 fHitVsSimChgVec.at(plane) = dir.make<TH2F>((
"HitVSimQ" +
std::to_string(plane)).c_str(),
"# e vs Hit SumADC;SumADC;# e", 50, 0., 1000., 50, 0., 100000.);
298 fHitVsSimIntVec.at(plane) = dir.make<TH2F>((
"HitVSimI" +
std::to_string(plane)).c_str(),
"# e vs Hit Integral;Integral;# e", 50, 0., 1000., 50, 0., 100000.);
299 fToteVHitEIntVec.at(plane) = dir.make<TH2F>((
"ToteVHite" +
std::to_string(plane)).c_str(),
"Tot e vs Hit e;Total #e;Hit #e", 50, 0., 100000., 50, 0., 100000.);
301 fWireEfficVec.at(plane) = dir.make<TProfile> ((
"WireEffic" +
std::to_string(plane)).c_str(),
"Wire Efficiency;# electrons", 50, 0., 100000., 0., 1.);
302 fWireEfficPHVec.at(plane) = dir.make<TProfile> ((
"WireEfficPH" +
std::to_string(plane)).c_str(),
"Wire Efficiency;# electrons", 50, 0., 20000., 0., 1.);
304 fHitEfficVec.at(plane) = dir.make<TProfile> ((
"HitEffic" +
std::to_string(plane)).c_str(),
"Hit Efficiency;Total # e-", 50, 0., 100000., 0., 1.);
305 fHitEfficPHVec.at(plane) = dir.make<TProfile> ((
"HitEfficPH" +
std::to_string(plane)).c_str(),
"Hit Efficiency;Max # e-", 50, 0., 20000., 0., 1.);
306 fHitEfficXZVec.at(plane) = dir.make<TProfile> ((
"HitEfficXZ" +
std::to_string(plane)).c_str(),
"Hit Efficiency;cos(thetaXZ)", 50, -1., 1., 0., 1.);
307 fHitEfficRMSVec.at(plane) = dir.make<TProfile> ((
"HitEfficRMS" +
std::to_string(plane)).c_str(),
"Hit Efficiency;MC Length (ticks)", 50, 0., 100., 0., 1.);
308 fCosXZvRMSVec.at(plane) = dir.make<TProfile> ((
"CosXZVRMS" +
std::to_string(plane)).c_str(),
"CosXZ v. RMS;Cos(XZ);Ticks", 50, -1., 1., 0., 100.);
309 fHitENEvXZVec.at(plane) = dir.make<TProfile2D>((
"HitENEvXZ" +
std::to_string(plane)).c_str(),
"XZ v # e;cos(thetaXZ);#electrons", 50, -1., 1.,
310 50, 0., 100000., 0., 1.);
316 void TrackHitEfficiencyAnalysis::initializeTuple(TTree* tree)
320 fTree->Branch(
"CryostataVec",
"std::vector<int>", &
fCryoVec);
355 void TrackHitEfficiencyAnalysis::clear()
const
390 void TrackHitEfficiencyAnalysis::fillHistograms(
const art::Event& event)
const
392 std::cout <<
" filling histos " << std::endl;
400 art::Handle< std::vector<sim::SimChannel>> simChannelHandle;
403 art::Handle< std::vector<simb::MCParticle>> mcParticleHandle;
407 if (!simChannelHandle.isValid() || simChannelHandle->empty() || !mcParticleHandle.isValid())
409 std::cout<<
"No sim channel information"<<std::endl;
412 std::cout<<
"There is sim channel information"<<std::endl;
418 using TDCToIDEMap = std::map<unsigned short, sim::IDE>;
419 using ChanToTDCToIDEMap = std::map<raw::ChannelID_t, TDCToIDEMap>;
420 using PartToChanToTDCToIDEMap = std::unordered_map<int, ChanToTDCToIDEMap>;
422 PartToChanToTDCToIDEMap partToChanToTDCToIDEMap;
425 for(
const auto& simChannel : *simChannelHandle)
427 for(
const auto& tdcide : simChannel.TDCIDEMap())
429 for(
const auto& ide : tdcide.second) partToChanToTDCToIDEMap[ide.trackID][simChannel.Channel()][tdcide.first] = ide;
436 using ChanToWireMap = std::unordered_map<raw::ChannelID_t,const recob::Wire*>;
438 ChanToWireMap channelToWireMap;
442 using ChanToRawDigitMap = std::unordered_map<raw::ChannelID_t,const raw::RawDigit*>;
444 ChanToRawDigitMap chanToRawDigitMap;
447 art::Handle< std::vector<int>> badChannelHandle;
456 art::Handle< std::vector<raw::RawDigit> > rawDigitHandle;
460 art::Handle< std::vector<recob::Wire> > wireHandle;;
465 if (!rawDigitHandle.isValid() || !wireHandle.isValid())
470 for(
const auto& wire : *wireHandle) channelToWireMap[wire.Channel()] = &wire;
472 for(
const auto& rawDigit : *rawDigitHandle) chanToRawDigitMap[rawDigit.Channel()] = &rawDigit;
477 using ChanToHitVecMap = std::unordered_map<raw::ChannelID_t,std::vector<const recob::Hit*>>;
479 ChanToHitVecMap channelToHitVec;
484 art::Handle< std::vector<recob::Hit> > hitHandle;
485 event.getByLabel(hitLabel, hitHandle);
487 for(
const auto&
hit : *hitHandle) channelToHitVec[
hit.Channel()].push_back(&
hit);
491 using TrackIDToMCParticleMap = std::unordered_map<int, const simb::MCParticle*>;
493 TrackIDToMCParticleMap trackIDToMCParticleMap;
495 for(
const auto& mcParticle : *mcParticleHandle)
496 trackIDToMCParticleMap[mcParticle.TrackId()] = &mcParticle;
499 auto const clockData = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(event);
501 std::vector<int> nSimChannelHitVec = {0,0,0};
502 std::vector<int> nRecobHitVec = {0,0,0};
503 std::vector<int> nFakeHitVec = {0,0,0};
504 std::vector<int> nSimulatedWiresVec = {0,0,0};
506 unsigned int lastwire=-1;
509 for(
const auto& partToChanInfo : partToChanToTDCToIDEMap)
512 TrackIDToMCParticleMap::const_iterator trackIDToMCPartItr = trackIDToMCParticleMap.find(partToChanInfo.first);
514 if (trackIDToMCPartItr == trackIDToMCParticleMap.end())
continue;
516 int trackPDGCode = trackIDToMCPartItr->second->PdgCode();
517 std::string processName = trackIDToMCPartItr->second->Process();
520 if (fabs(trackPDGCode) != 13 || processName !=
"primary")
continue;
523 Eigen::Vector3f partStartPos(trackIDToMCPartItr->second->Vx(),trackIDToMCPartItr->second->Vy(),trackIDToMCPartItr->second->Vz());
524 Eigen::Vector3f partStartDir(trackIDToMCPartItr->second->Px(),trackIDToMCPartItr->second->Py(),trackIDToMCPartItr->second->Pz());
526 partStartDir.normalize();
528 Eigen::Vector2f partStartDirVecXZ(partStartDir[0],partStartDir[2]);
530 partStartDirVecXZ.normalize();
534 std::vector<Eigen::Vector3f> lastPositionVec = {partStartPos,partStartPos,partStartPos};
536 for(
const auto& chanToTDCToIDEMap : partToChanInfo.second)
545 std::cout <<
"*** skipping bad channel with status: " << chanFilt.
Status(chanToTDCToIDEMap.first) <<
" for channel: " << chanToTDCToIDEMap.first <<
", plane: " << wids[0].Plane <<
", wire: " << wids[0].Wire << std::endl;
552 if (badChannelHandle.isValid())
555 std::vector<int>::const_iterator badItr = std::find(badChannelHandle->begin(),badChannelHandle->end(),chanToTDCToIDEMap.first);
557 if (badItr != badChannelHandle->end())
continue;
578 TDCToIDEMap tdcToIDEMap = chanToTDCToIDEMap.second;
579 float totalElectrons(0.);
580 float maxElectrons(0.);
581 unsigned short maxElectronsTDC(0);
582 int nMatchedWires(0);
590 unsigned int plane = wids[0].Plane;
591 unsigned int wire = wids[0].Wire;
593 Eigen::Vector3f avePosition(0.,0.,0.);
595 if(wire!=lastwire) nSimulatedWiresVec[plane]++;
598 for(
const auto& ideVal : tdcToIDEMap)
600 totalElectrons += ideVal.second.numElectrons;
602 if (maxElectrons < ideVal.second.numElectrons)
604 maxElectrons = ideVal.second.numElectrons;
605 maxElectronsTDC = ideVal.first;
608 avePosition += Eigen::Vector3f(ideVal.second.x,ideVal.second.y,ideVal.second.z);
613 avePosition /= float(tdcToIDEMap.size());
615 Eigen::Vector3f partDirVec = avePosition - lastPositionVec[plane];
617 partDirVec.normalize();
619 lastPositionVec[plane] = avePosition;
622 Eigen::Vector2f projPairDirVec(partDirVec[0],partDirVec[2]);
624 projPairDirVec.normalize();
626 float cosThetaXZ = projPairDirVec[0];
628 totalElectrons = std::min(totalElectrons,
float(99900.));
633 nSimChannelHitVec[plane]++;
635 unsigned short startTDC = tdcToIDEMap.begin()->first;
636 unsigned short stopTDC = tdcToIDEMap.rbegin()->first;
639 unsigned short startTick = clockData.TPCTDC2Tick(startTDC) +
fOffsetVec[plane];
640 unsigned short stopTick = clockData.TPCTDC2Tick(stopTDC) +
fOffsetVec[plane];
641 unsigned short maxETick = clockData.TPCTDC2Tick(maxElectronsTDC) +
fOffsetVec[plane];
646 float nElectronsTotalBest(0.);
647 float hitSummedADCBest(0.);
648 float hitIntegralBest(0.);
649 float hitPeakTimeBest(0.);
650 float hitPeakAmpBest(-100.);
651 float hitRMSBest(0.);
652 int hitMultiplicityBest(0);
653 int hitLocalIndexBest(0);
654 float hitGoodnessBest(0.);
655 int hitNumDegreesBest(0);
656 float hitBaselineBest(0.);
657 float hitSnippetLenBest(0.);
658 unsigned short hitStopTickBest(0);
659 unsigned short hitStartTickBest(0);
662 ChanToWireMap::const_iterator wireItr = channelToWireMap.find(chanToTDCToIDEMap.first);
664 if (wireItr != channelToWireMap.end())
670 for(
const auto& range : signalROI.
get_ranges())
677 const std::vector<float>& signal = range.data();
682 if (roiFirstBinTick > stopTick || roiLastBinTick < startTick)
continue;
684 wireRangePtr = ⦥
699 ChanToHitVecMap::iterator hitIter = channelToHitVec.find(chanToTDCToIDEMap.first);
701 if (hitIter != channelToHitVec.end())
707 for(
const auto&
hit : hitIter->second)
709 unsigned short hitStartTick =
hit->PeakTime() -
fSigmaVec[plane] *
hit->RMS();
710 unsigned short hitStopTick =
hit->PeakTime() +
fSigmaVec[plane] *
hit->RMS();
713 if (hitStartTick > stopTick || hitStopTick < startTick)
715 nFakeHitVec[plane]++;
720 float hitHeight =
hit->PeakAmplitude();
723 if (hitHeight < hitPeakAmpBest)
continue;
725 hitPeakAmpBest = hitHeight;
727 hitStartTickBest = hitStartTick;
728 hitStopTickBest = hitStopTick;
735 nElectronsTotalBest = 0.;
736 hitPeakTimeBest = bestHit->
PeakTime();
737 hitIntegralBest = bestHit->
Integral();
739 hitRMSBest = bestHit->
RMS();
745 hitBaselineBest = 0.;
750 for(
unsigned short tick = hitStartTickBest;
tick <= hitStopTickBest;
tick++)
752 unsigned short hitTDC = clockData.TPCTick2TDC(
tick -
fOffsetVec[plane]);
754 TDCToIDEMap::iterator ideIterator = tdcToIDEMap.find(hitTDC);
756 if (ideIterator != tdcToIDEMap.end()) nElectronsTotalBest += ideIterator->second.numElectrons;
760 if (nMatchedHits > 0)
762 float chgRatioADC = hitSummedADCBest > 0. ? totalElectrons / hitSummedADCBest : 0.;
763 float chgRatioInt = hitIntegralBest > 0. ? totalElectrons / hitIntegralBest : 0.;
769 fHitVsSimChgVec[plane]->Fill(std::min(hitSummedADCBest,
float(999.)), totalElectrons, 1.);
770 fHitVsSimIntVec[plane]->Fill(std::min(hitIntegralBest,
float(999.)), totalElectrons, 1.);
771 fToteVHitEIntVec[plane]->Fill(std::min(totalElectrons,
float(99999.)),std::min(nElectronsTotalBest,
float(99999.)),1.);
775 fHitNumTDCVec[plane]->Fill(std::min(
float(hitStopTickBest - hitStartTickBest),
float(99.5)), 1.);
776 fSnippetLenVec[plane]->Fill(std::min(hitSnippetLenBest,
float(99.5)), 1.);
779 nRecobHitVec[plane]++;
782 else if (rejectedHit)
787 mf::LogDebug(
"TrackHitEfficiencyAnalysis") <<
"**> TPC: " << rejectedHit->
WireID().
TPC <<
", Plane " << rejectedHit->
WireID().
Plane <<
", wire: " << rejectedHit->
WireID().
Wire <<
", hit startstop tick: " << hitStartTick <<
"/" << hitStopTick <<
", start/stop ticks: " << startTick <<
"/" << stopTick << std::endl;
788 mf::LogDebug(
"TrackHitEfficiencyAnalysis") <<
" TPC/Plane/Wire: " << wids[0].TPC <<
"/" << plane <<
"/" << wids[0].Wire <<
", Track # hits: " << partToChanInfo.second.size() <<
", # hits: "<< hitIter->second.size() <<
", # electrons: " << totalElectrons <<
", pulse Height: " << rejectedHit->
PeakAmplitude() <<
", charge: " << rejectedHit->
Integral() <<
", " <<rejectedHit->
SummedADC() << std::endl;
792 mf::LogDebug(
"TrackHitEfficiencyAnalysis") <<
"==> No match, TPC/Plane/Wire: " <<
"/" << wids[0].TPC <<
"/" << wids[0].Plane <<
"/" << wids[0].Wire <<
", # electrons: " << totalElectrons <<
", startTick: " << startTick <<
", stopTick: " << stopTick << std::endl;
798 fWireEfficVec.at(plane)->Fill(totalElectrons, std::min(nMatchedWires,1), 1.);
799 fWireEfficPHVec.at(plane)->Fill(maxElectrons, std::min(nMatchedWires,1), 1.);
801 float matchHit = std::min(nMatchedHits,1);
802 float snippetLen = std::min(
float(stopTick - startTick),
float(99.5));
805 fHitEfficVec[plane]->Fill(totalElectrons, matchHit, 1.);
811 fHitENEvXZVec[plane]->Fill(cosThetaXZ, totalElectrons, matchHit);
815 fCryoVec.push_back(wids[0].Cryostat);
848 if (nSimChannelHitVec[idx] > 10)
850 float hitEfficiency = float(nRecobHitVec[idx]) / float(nSimChannelHitVec[idx]);
853 fNRecobHitVec[idx]->Fill(std::min(nRecobHitVec[idx],1999), 1.);
854 fNFakeHitVec[idx]->Fill(nFakeHitVec[idx]/(
float)nSimulatedWiresVec[idx],1.);
863 void TrackHitEfficiencyAnalysis::endJob(
int numEvents)
868 DEFINE_ART_CLASS_TOOL(TrackHitEfficiencyAnalysis)
short int LocalIndex() const
How well do we believe we know this hit?
std::vector< float > fHitIntegralVec
std::vector< TH1F * > fNFakeHitVec
std::vector< TProfile * > fHitEfficXZVec
std::vector< TH1F * > fTotalElectronsHistVec
std::vector< float > fSigmaVec
Window size for matching to SimChannels.
std::vector< TH1F * > fSimDivHitChgVec
std::vector< int > fStopTickVec
Utilities related to art service access.
BEGIN_PROLOG true icarus_rawdigitfilter FilterTools FilterPlane1 Plane
std::vector< TProfile * > fHitEfficPHVec
std::vector< float > fHitPeakAmpVec
std::vector< float > fHitBaselinevec
std::vector< int > fPlaneVec
geo::WireID WireID() const
float RMS() const
RMS of the hit shape, in tick units.
Declaration of signal hit object.
std::vector< TH1F * > fHitEfficiencyVec
std::vector< int > fNumDegreesVec
std::vector< TH1F * > fSnippetLenVec
std::vector< TH2F * > fToteVHitEIntVec
std::vector< int > fMaxETickVec
int DegreesOfFreedom() const
std::vector< geo::WireID > ChannelToWire(raw::ChannelID_t const channel) const
Returns a list of wires connected to the specified TPC channel.
Definition of basic raw digits.
void configure(fhicl::ParameterSet const &pset) override
const range_list_t & get_ranges() const
Returns the internal list of non-void ranges.
std::vector< int > fStartTickVec
float Integral() const
Integral under the calibrated signal waveform of the hit, in tick x ADC units.
std::vector< TProfile * > fWireEfficVec
WireID_t Wire
Index of the wire within its plane.
std::vector< float > fMaxElectronsVec
std::vector< TH2F * > fHitVsSimIntVec
float GoodnessOfFit() const
Degrees of freedom in the determination of the hit signal shape (-1 by default)
std::vector< TH1F * > fSimNumTDCVec
short int Multiplicity() const
How many hits could this one be shared with.
art::InputTag fBadChannelProducerLabel
std::vector< int > fOffsetVec
Allow offsets for each plane.
float PeakAmplitude() const
The estimated amplitude of the hit at its peak, in ADC units.
void initializeHists(art::ServiceHandle< art::TFileService > &, const std::string &) override
Interface for initializing the histograms to be filled.
int TDCtick_t
Type representing a TDC tick.
std::vector< TH1F * > fDeltaMidTDCVec
void endJob(int numEvents) override
Interface for method to executve at the end of run processing.
std::vector< TH1F * > fHitIntegralHistVec
std::vector< art::Ptr< recob::Hit >> HitPtrVec
std::vector< TProfile * > fHitEfficRMSVec
std::vector< TH1F * > fHitNumTDCVec
art::InputTag fMCParticleProducerLabel
unsigned int Nplanes(unsigned int tpc=0, unsigned int cstat=0) const
Returns the total number of wire planes in the specified TPC.
std::vector< int > fWireVec
void fillHistograms(const art::Event &) const override
Interface for filling histograms.
std::vector< TProfile * > fWireEfficPHVec
void initializeTuple(TTree *) override
Interface for initializing the tuple variables.
std::vector< int > fTPCVec
std::vector< TH1F * > fHitElectronsVec
std::vector< int > fHitMultiplicityVec
std::vector< float > fHitGoodnessVec
std::vector< TH1F * > fNMatchedHitVec
tick_as<> tick
Tick number, represented by std::ptrdiff_t.
virtual Status_t Status(raw::ChannelID_t channel) const
Returns a status integer with arbitrary meaning.
std::vector< TH1F * > fHitPulseHeightVec
art::InputTag fSimChannelProducerLabel
std::string fLocalDirName
Fraction for truncated mean.
Class providing information about the quality of channels.
std::vector< TProfile * > fHitEfficVec
PlaneID_t Plane
Index of the plane within its TPC.
std::vector< int > fNMatchedWires
Description of geometry of one entire detector.
std::map< size_t, HitPtrVec > ViewHitMap
std::vector< art::InputTag > fHitProducerLabelVec
raw::TDCtick_t StartTick() const
Initial tdc tick for hit.
~TrackHitEfficiencyAnalysis()
Destructor.
std::vector< float > fHitPeakTimeVec
std::vector< TProfile * > fCosXZvRMSVec
raw::TDCtick_t EndTick() const
Final tdc tick for hit.
std::vector< art::InputTag > fRawDigitProducerLabelVec
std::vector< TH1F * > fHitPulseWidthVec
std::vector< int > fHitLocalIndexVec
std::vector< TH1F * > fHitSumADCVec
std::vector< int > fCryoVec
float PeakTime() const
Time of the signal peak, in tick units.
TrackHitEfficiencyAnalysis(fhicl::ParameterSet const &pset)
Constructor.
std::vector< art::InputTag > fWireProducerLabelVec
std::vector< TProfile2D * > fHitENEvXZVec
std::vector< TH1F * > fNSimChannelHitsVec
std::string to_string(WindowPattern const &pattern)
std::vector< TH2F * > fHitVsSimChgVec
std::vector< int > fHitStopTickVec
object containing MC truth information necessary for making RawDigits and doing back tracking ...
std::vector< int > fHitStartTickVec
Interface for experiment-specific channel quality info provider.
std::vector< TH1F * > fNRecobHitVec
int fMinAllowedChanStatus
Don't consider channels with lower status.
std::vector< float > fPartDirY
Range class, with range and data.
float SummedADC() const
The sum of calibrated ADC counts of the hit (0. by default)
Declaration of basic channel signal object.
std::vector< TH1F * > fSimDivHitChg1Vec
2D representation of charge deposited in the TDC/wire plane
art::ServiceHandle< art::TFileService > tfs
TPCID_t TPC
Index of the TPC within its cryostat.
Interface for experiment-specific service for channel quality info.
std::map< int, ViewHitMap > TrackViewHitMap
std::vector< float > fPartDirZ
std::vector< float > fHitSummedADCVec
std::vector< float > fTotalElectronsVec
std::vector< float > fHitPeakRMSVec
std::vector< float > fPartDirX
art framework interface to geometry description
BEGIN_PROLOG could also be cout
std::vector< TH1F * > fMaxElectronsHistVec
std::vector< int > fNMatchedHits
const geo::GeometryCore * fGeometry
pointer to Geometry service