All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
icaruscode/icaruscode/Analysis/tools/TrackHitEfficiencyAnalysis_tool.cc
Go to the documentation of this file.
2 
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"
10 
12 #include "larcore/CoreUtils/ServiceUtil.h" // lar::providerFrom()
16 
21 #include "nusimdata/SimulationBase/MCParticle.h"
22 
23 // Eigen
24 #include <Eigen/Dense>
25 
26 #include "TH1.h"
27 #include "TH2.h"
28 #include "TProfile.h"
29 #include "TProfile2D.h"
30 #include "TTree.h"
31 
32 #include <cmath>
33 #include <algorithm>
34 
36 {
37  ////////////////////////////////////////////////////////////////////////
38  //
39  // Class: TrackHitEfficiencyAnalysis
40  // Module Type: producer
41  // File: TrackHitEfficiencyAnalysis.cc
42  //
43  // The intent of this module is to provide methods for
44  // "analyzing" hits on waveforms
45  //
46  // Configuration parameters:
47  //
48  // TruncMeanFraction - the fraction of waveform bins to discard when
49  //
50  // Created by Tracy Usher (usher@slac.stanford.edu) on February 19, 2016
51  //
52  ////////////////////////////////////////////////////////////////////////
53 
54 // The following typedefs will, obviously, be useful
55 using HitPtrVec = std::vector<art::Ptr<recob::Hit>>;
56 using ViewHitMap = std::map<size_t,HitPtrVec>;
57 using TrackViewHitMap = std::map<int,ViewHitMap>;
58 
60 {
61 public:
62  /**
63  * @brief Constructor
64  *
65  * @param pset
66  */
67  explicit TrackHitEfficiencyAnalysis(fhicl::ParameterSet const & pset);
68 
69  /**
70  * @brief Destructor
71  */
73 
74  // provide for initialization
75  void configure(fhicl::ParameterSet const & pset) override;
76 
77  /**
78  * @brief Interface for initializing the histograms to be filled
79  *
80  * @param TFileService handle to the TFile service
81  * @param string subdirectory to store the hists in
82  */
83  void initializeHists(art::ServiceHandle<art::TFileService>&, const std::string&) override;
84 
85  /**
86  * @brief Interface for initializing the tuple variables
87  *
88  * @param TTree pointer to a TTree object to which to add variables
89  */
90  void initializeTuple(TTree*) override;
91 
92  /**
93  * @brief Interface for method to executve at the end of run processing
94  *
95  * @param int number of events to use for normalization
96  */
97  void endJob(int numEvents) override;
98 
99  /**
100  * @brief Interface for filling histograms
101  */
102  void fillHistograms(const art::Event&) const override;
103 
104 private:
105 
106  // Clear mutable variables
107  void clear() const;
108 
109  // Fcl parameters.
110  std::vector<art::InputTag> fRawDigitProducerLabelVec;
111  std::vector<art::InputTag> fWireProducerLabelVec;
112  std::vector<art::InputTag> fHitProducerLabelVec;
117  std::string fLocalDirName; ///< Fraction for truncated mean
118  std::vector<int> fOffsetVec; ///< Allow offsets for each plane
119  std::vector<float> fSigmaVec; ///< Window size for matching to SimChannels
120  int fMinAllowedChanStatus; ///< Don't consider channels with lower status
122 
123  // Pointers to the histograms we'll create.
124  std::vector<TH1F*> fTotalElectronsHistVec;
125  std::vector<TH1F*> fMaxElectronsHistVec;
126  std::vector<TH1F*> fHitElectronsVec;
127  std::vector<TH1F*> fHitSumADCVec;
128  std::vector<TH1F*> fHitIntegralHistVec;
129  std::vector<TH1F*> fHitPulseHeightVec;
130  std::vector<TH1F*> fHitPulseWidthVec;
131  std::vector<TH1F*> fSimNumTDCVec;
132  std::vector<TH1F*> fHitNumTDCVec;
133  std::vector<TH1F*> fSnippetLenVec;
134  std::vector<TH1F*> fNMatchedHitVec;
135  std::vector<TH1F*> fDeltaMidTDCVec;
136  std::vector<TProfile*> fWireEfficVec;
137  std::vector<TProfile*> fWireEfficPHVec;
138  std::vector<TProfile*> fHitEfficVec;
139  std::vector<TProfile*> fHitEfficPHVec;
140  std::vector<TProfile*> fHitEfficXZVec;
141  std::vector<TProfile*> fHitEfficRMSVec;
142  std::vector<TProfile*> fCosXZvRMSVec;
143  std::vector<TProfile2D*> fHitENEvXZVec;
144  std::vector<TH1F*> fSimDivHitChgVec;
145  std::vector<TH1F*> fSimDivHitChg1Vec;
146  std::vector<TH2F*> fHitVsSimChgVec;
147  std::vector<TH2F*> fHitVsSimIntVec;
148  std::vector<TH2F*> fToteVHitEIntVec;
149 
150  std::vector<TH1F*> fNSimChannelHitsVec;
151  std::vector<TH1F*> fNRecobHitVec;
152  //std::vector<TH1F*> fNRejectedHitVec;
153  std::vector<TH1F*> fHitEfficiencyVec;
154 
155  std::vector<TH1F*> fNFakeHitVec;
156 
157  // TTree variables
158  mutable TTree* fTree;
159 
160  mutable std::vector<int> fTPCVec;
161  mutable std::vector<int> fCryoVec;
162  mutable std::vector<int> fPlaneVec;
163  mutable std::vector<int> fWireVec;
164 
165  mutable std::vector<float> fTotalElectronsVec;
166  mutable std::vector<float> fMaxElectronsVec;
167  mutable std::vector<int> fStartTickVec;
168  mutable std::vector<int> fStopTickVec;
169  mutable std::vector<int> fMaxETickVec;
170  mutable std::vector<float> fPartDirX;
171  mutable std::vector<float> fPartDirY;
172  mutable std::vector<float> fPartDirZ;
173 
174  mutable std::vector<int> fNMatchedWires;
175  mutable std::vector<int> fNMatchedHits;
176 
177  mutable std::vector<float> fHitPeakTimeVec;
178  mutable std::vector<float> fHitPeakAmpVec;
179  mutable std::vector<float> fHitPeakRMSVec;
180  mutable std::vector<float> fHitBaselinevec;
181  mutable std::vector<float> fHitSummedADCVec;
182  mutable std::vector<float> fHitIntegralVec;
183  mutable std::vector<int> fHitStartTickVec;
184  mutable std::vector<int> fHitStopTickVec;
185  mutable std::vector<int> fHitMultiplicityVec;
186  mutable std::vector<int> fHitLocalIndexVec;
187  mutable std::vector<float> fHitGoodnessVec;
188  mutable std::vector<int> fNumDegreesVec;
189 
190  // Useful services, keep copies for now (we can update during begin run periods)
191  const geo::GeometryCore* fGeometry; ///< pointer to Geometry service
192 };
193 
194 //----------------------------------------------------------------------------
195 /// Constructor.
196 ///
197 /// Arguments:
198 ///
199 /// pset - Fcl parameters.
200 ///
201 TrackHitEfficiencyAnalysis::TrackHitEfficiencyAnalysis(fhicl::ParameterSet const & pset) : fTree(nullptr)
202 {
203  fGeometry = lar::providerFrom<geo::Geometry>();
204 
205  configure(pset);
206 
207  // Report.
208  mf::LogInfo("TrackHitEfficiencyAnalysis") << "TrackHitEfficiencyAnalysis configured\n";
209 }
210 
211 //----------------------------------------------------------------------------
212 /// Destructor.
213 TrackHitEfficiencyAnalysis::~TrackHitEfficiencyAnalysis()
214 {}
215 
216 //----------------------------------------------------------------------------
217 /// Reconfigure method.
218 ///
219 /// Arguments:
220 ///
221 /// pset - Fcl parameter set.
222 ///
223 void TrackHitEfficiencyAnalysis::configure(fhicl::ParameterSet const & pset)
224 {
225  fRawDigitProducerLabelVec = pset.get< std::vector<art::InputTag>>("RawDigitLabelVec", std::vector<art::InputTag>() = {"rawdigitfilter"});
226  fWireProducerLabelVec = pset.get< std::vector<art::InputTag>>("WireModuleLabelVec", std::vector<art::InputTag>() = {"decon1droi"});
227  fHitProducerLabelVec = pset.get< std::vector<art::InputTag>>("HitModuleLabelVec", std::vector<art::InputTag>() = {"gauss"});
228  fMCParticleProducerLabel = pset.get< art::InputTag >("MCParticleLabel", "largeant");
229  fSimChannelProducerLabel = pset.get< art::InputTag >("SimChannelLabel", "largeant");
230  fBadChannelProducerLabel = pset.get< art::InputTag >("BadChannelLabel", "simnfspl1:badchannels");
231  fUseBadChannelDB = pset.get< bool >("UseBadChannelDB", true);
232  fLocalDirName = pset.get<std::string >("LocalDirName", std::string("wow"));
233  fOffsetVec = pset.get<std::vector<int> >("OffsetVec", std::vector<int>()={0,0,0});
234  fSigmaVec = pset.get<std::vector<float> >("SigmaVec", std::vector<float>()={1.,1.,1.});
235  fMinAllowedChanStatus = pset.get< int >("MinAllowedChannelStatus");
236  fSimChannelMinEnergy = pset.get<float >("SimChannelMinEnergy", std::numeric_limits<float>::epsilon());
237 }
238 
239 //----------------------------------------------------------------------------
240 /// Begin job method.
241 void TrackHitEfficiencyAnalysis::initializeHists(art::ServiceHandle<art::TFileService>& tfs, const std::string& dirName)
242 {
243  // Make a directory for these histograms
244  art::TFileDirectory dir = tfs->mkdir(dirName.c_str());
245 
249  fHitSumADCVec.resize(fGeometry->Nplanes());
253  fSimNumTDCVec.resize(fGeometry->Nplanes());
254  fHitNumTDCVec.resize(fGeometry->Nplanes());
255  fSnippetLenVec.resize(fGeometry->Nplanes());
256  fNMatchedHitVec.resize(fGeometry->Nplanes());
257  fDeltaMidTDCVec.resize(fGeometry->Nplanes());
258  fHitVsSimChgVec.resize(fGeometry->Nplanes());
259  fHitVsSimIntVec.resize(fGeometry->Nplanes());
260  fCosXZvRMSVec.resize(fGeometry->Nplanes());
262  fHitENEvXZVec.resize(fGeometry->Nplanes());
264  fNRecobHitVec.resize(fGeometry->Nplanes());
265  fNFakeHitVec.resize(fGeometry->Nplanes());
269 
270  fWireEfficVec.resize(fGeometry->Nplanes());
271  fWireEfficPHVec.resize(fGeometry->Nplanes());
272 
273  fHitEfficVec.resize(fGeometry->Nplanes());
274  fHitEfficPHVec.resize(fGeometry->Nplanes());
275  fHitEfficXZVec.resize(fGeometry->Nplanes());
276  fHitEfficRMSVec.resize(fGeometry->Nplanes());
277 
278  for(size_t plane = 0; plane < fGeometry->Nplanes(); plane++)
279  {
280  fTotalElectronsHistVec.at(plane) = dir.make<TH1F>(("TotalElecs" + std::to_string(plane)).c_str(), ";Total # electrons", 250, 0., 100000.);
281  fMaxElectronsHistVec.at(plane) = dir.make<TH1F>(("MaxElecs" + std::to_string(plane)).c_str(), ";Max # electrons", 250, 0., 20000.);
282  fHitElectronsVec.at(plane) = dir.make<TH1F>(("HitElecs" + std::to_string(plane)).c_str(), ";# e- in Hit Range", 250, 0., 100000.);
283  fHitSumADCVec.at(plane) = dir.make<TH1F>(("SumADC" + std::to_string(plane)).c_str(), "Hit Sum ADC", 200, 0., 1000.);
284  fHitIntegralHistVec.at(plane) = dir.make<TH1F>(("Integral" + std::to_string(plane)).c_str(), "Hit Integral", 200, 0., 1000.);
285  fHitPulseHeightVec.at(plane) = dir.make<TH1F>(("PulseHeight" + std::to_string(plane)).c_str(), "Hit PH (ADC)", 200, 0., 150.);
286  fHitPulseWidthVec.at(plane) = dir.make<TH1F>(("PulseWidth" + std::to_string(plane)).c_str(), "Hit PulseWidth;Hit RMS", 100, 0., 20.);
287  fSimNumTDCVec.at(plane) = dir.make<TH1F>(("SimNumTDC" + std::to_string(plane)).c_str(), ";TDC ticks", 100, 0., 100.);
288  fHitNumTDCVec.at(plane) = dir.make<TH1F>(("HitNumTDC" + std::to_string(plane)).c_str(), ";ticks", 100, 0., 100.);
289  fSnippetLenVec.at(plane) = dir.make<TH1F>(("SnippetLen" + std::to_string(plane)).c_str(), ";ticks", 100, 0., 100.);
290  fNMatchedHitVec.at(plane) = dir.make<TH1F>(("NMatched" + std::to_string(plane)).c_str(), ";# hits", 20, 0., 20.);
291  fDeltaMidTDCVec.at(plane) = dir.make<TH1F>(("DeltaMid" + std::to_string(plane)).c_str(), ";# hits", 50, -25., 25.);
292  fNSimChannelHitsVec.at(plane) = dir.make<TH1F>(("NSimChan" + std::to_string(plane)).c_str(), ";# hits", 100, 0., 1200.);
293  fNRecobHitVec.at(plane) = dir.make<TH1F>(("NRecobHit" + std::to_string(plane)).c_str(), ";# hits", 100, 0., 1200.);
294  fNFakeHitVec.at(plane) = dir.make<TH1F>(("NFakeHit" + std::to_string(plane)).c_str(), ";# hits", 100, 0., 50.);
295  fHitEfficiencyVec.at(plane) = dir.make<TH1F>(("PlnEffic" + std::to_string(plane)).c_str(), ";# hits", 101, 0., 1.01);
296  fSimDivHitChgVec.at(plane) = dir.make<TH1F>(("SimDivHit" + std::to_string(plane)).c_str(), ";# e / SummedADC", 200, 0., 200.);
297  fSimDivHitChg1Vec.at(plane) = dir.make<TH1F>(("SimDivHit1" + std::to_string(plane)).c_str(), ";# e / Integral", 200, 0., 200.);
298 
299  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.);
300  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.);
301  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.);
302 
303  fWireEfficVec.at(plane) = dir.make<TProfile> (("WireEffic" + std::to_string(plane)).c_str(), "Wire Efficiency;# electrons", 50, 0., 100000., 0., 1.);
304  fWireEfficPHVec.at(plane) = dir.make<TProfile> (("WireEfficPH" + std::to_string(plane)).c_str(), "Wire Efficiency;# electrons", 50, 0., 20000., 0., 1.);
305 
306  fHitEfficVec.at(plane) = dir.make<TProfile> (("HitEffic" + std::to_string(plane)).c_str(), "Hit Efficiency;Total # e-", 50, 0., 100000., 0., 1.);
307  fHitEfficPHVec.at(plane) = dir.make<TProfile> (("HitEfficPH" + std::to_string(plane)).c_str(), "Hit Efficiency;Max # e-", 50, 0., 20000., 0., 1.);
308  fHitEfficXZVec.at(plane) = dir.make<TProfile> (("HitEfficXZ" + std::to_string(plane)).c_str(), "Hit Efficiency;cos(thetaXZ)", 50, -1., 1., 0., 1.);
309  fHitEfficRMSVec.at(plane) = dir.make<TProfile> (("HitEfficRMS" + std::to_string(plane)).c_str(), "Hit Efficiency;MC Length (ticks)", 50, 0., 100., 0., 1.);
310  fCosXZvRMSVec.at(plane) = dir.make<TProfile> (("CosXZVRMS" + std::to_string(plane)).c_str(), "CosXZ v. RMS;Cos(XZ);Ticks", 50, -1., 1., 0., 100.);
311  fHitENEvXZVec.at(plane) = dir.make<TProfile2D>(("HitENEvXZ" + std::to_string(plane)).c_str(), "XZ v # e;cos(thetaXZ);#electrons", 50, -1., 1.,
312  50, 0., 100000., 0., 1.);
313  }
314 
315  return;
316 }
317 
318 void TrackHitEfficiencyAnalysis::initializeTuple(TTree* tree)
319 {
320  fTree = tree;
321 
322  fTree->Branch("CryostataVec", "std::vector<int>", &fCryoVec);
323  fTree->Branch("TPCVec", "std::vector<int>", &fTPCVec);
324  fTree->Branch("PlaneVec", "std::vector<int>", &fPlaneVec);
325  fTree->Branch("WireVec", "std::vector<int>", &fWireVec);
326 
327  fTree->Branch("TotalElectronsVec", "std::vector<float>", &fTotalElectronsVec);
328  fTree->Branch("MaxElectronsVec", "std::vector<float>", &fMaxElectronsVec);
329  fTree->Branch("StartTick", "std::vector<int>", &fStartTickVec);
330  fTree->Branch("StopTick", "std::vector<int>", &fStopTickVec);
331  fTree->Branch("MaxETick", "std::vector<int>", &fMaxETickVec);
332  fTree->Branch("PartDirX", "std::vector<float>", &fPartDirX);
333  fTree->Branch("PartDirY", "std::vector<float>", &fPartDirY);
334  fTree->Branch("PartDirZ", "std::vector<float>", &fPartDirZ);
335 
336  fTree->Branch("NMatchedWires", "std::vector<int>", &fNMatchedWires);
337  fTree->Branch("NMatchedHits", "std::vector<int>", &fNMatchedHits);
338 
339  fTree->Branch("HitPeakTimeVec", "std::vector<float>", &fHitPeakTimeVec);
340  fTree->Branch("HitPeakAmpVec", "std::vector<float>", &fHitPeakAmpVec);
341  fTree->Branch("HitPeakRMSVec", "std::vector<float>", &fHitPeakRMSVec);
342  fTree->Branch("HitBaselineVec", "std::vector<float>", &fHitBaselinevec);
343  fTree->Branch("HitSummedADCVec", "std::vector<float>", &fHitSummedADCVec);
344  fTree->Branch("HitIntegralVec", "std::vector<float>", &fHitIntegralVec);
345  fTree->Branch("HitStartTickVec", "std::vector<int>", &fHitStartTickVec);
346  fTree->Branch("HitStopTickVec", "std::vector<int>", &fHitStopTickVec);
347  fTree->Branch("HitMultiplicity", "std::vector<int>", &fHitMultiplicityVec);
348  fTree->Branch("HitLocalIndex", "std::vector<int>", &fHitLocalIndexVec);
349  fTree->Branch("HitGoodness", "std::vector<float>", &fHitGoodnessVec);
350  fTree->Branch("HitNumDegrees", "std::vector<int>", &fNumDegreesVec);
351 
352  clear();
353 
354  return;
355 }
356 
357 void TrackHitEfficiencyAnalysis::clear() const
358 {
359  fTPCVec.clear();
360  fCryoVec.clear();
361  fPlaneVec.clear();
362  fWireVec.clear();
363 
364  fTotalElectronsVec.clear();
365  fMaxElectronsVec.clear();
366  fStartTickVec.clear();
367  fStopTickVec.clear();
368  fMaxETickVec.clear();
369  fPartDirX.clear();
370  fPartDirY.clear();
371  fPartDirZ.clear();
372 
373  fNMatchedWires.clear();
374  fNMatchedHits.clear();
375 
376  fHitPeakTimeVec.clear();
377  fHitPeakAmpVec.clear();
378  fHitPeakRMSVec.clear();
379  fHitBaselinevec.clear();
380  fHitSummedADCVec.clear();
381  fHitIntegralVec.clear();
382  fHitStartTickVec.clear();
383  fHitStopTickVec.clear();
384  fHitMultiplicityVec.clear();
385  fHitLocalIndexVec.clear();
386  fHitGoodnessVec.clear();
387  fNumDegreesVec.clear();
388 
389  return;
390 }
391 
392 void TrackHitEfficiencyAnalysis::fillHistograms(const art::Event& event) const
393 {
394  // std::cout << " filling histos " << std::endl;
395  // Basic assumption is that the producer label vecs for RawDigits and Wire data are
396  // all the same length and in the same order. Here we just check for length
397  if (fRawDigitProducerLabelVec.size() != fWireProducerLabelVec.size()) return;
398 
399  // Always clear the tuple
400  clear();
401 
402  art::Handle< std::vector<sim::SimChannel>> simChannelHandle;
403  event.getByLabel(fSimChannelProducerLabel, simChannelHandle);
404 
405  art::Handle< std::vector<simb::MCParticle>> mcParticleHandle;
406  event.getByLabel(fMCParticleProducerLabel, mcParticleHandle);
407 
408  // If there is no sim channel informaton then exit
409  if (!simChannelHandle.isValid() || simChannelHandle->empty() || !mcParticleHandle.isValid()) return;
410 
411  // There are several things going on here... for each channel we have particles (track id's) depositing energy in a range to ticks
412  // So... for each channel we want to build a structure that relates particles to tdc ranges and deposited energy (or electrons)
413  // Here is a complicated structure:
414 // using TDCToIDEMap = std::map<unsigned short, sim::IDE>; // We need this one in order
415 // using ChanToTDCToIDEMap = std::map<raw::ChannelID_t, TDCToIDEMap>;
416  using TDCIDEPair = std::pair<unsigned short, const sim::IDE*>;
417  using TickTDCIDEVec = std::vector<TDCIDEPair>;
418  using ChanToTDCIDEMap = std::unordered_map<raw::ChannelID_t,TickTDCIDEVec>;
419  using PartToChanToTDCToIDEMap = std::unordered_map<int, ChanToTDCIDEMap>;
420 
421  PartToChanToTDCToIDEMap partToChanToTDCToIDEMap;
422 
423  // Build out the above data structure
424  for(const auto& simChannel : *simChannelHandle)
425  {
426  raw::ChannelID_t channel = simChannel.Channel();
427 
428  for(const auto& tdcide : simChannel.TDCIDEMap())
429  {
430  for(const auto& ide : tdcide.second) //chanToTDCToIDEMap[simChannel.Channel()][tdcide.first] = ide;
431  {
432  if (ide.energy < fSimChannelMinEnergy) continue;
433 
434  partToChanToTDCToIDEMap[ide.trackID][channel].emplace_back(tdcide.first,&ide);
435 
436  if (ide.energy < std::numeric_limits<float>::epsilon()) mf::LogDebug("SpacePointAnalysis") << ">> epsilon simchan deposited energy: " << ide.energy << std::endl;
437  }
438  }
439  }
440 
441  // what needs to be done?
442  // First we define a straightforward channel to Wire map so we can look up a given
443  // channel's Wire data as we loop over SimChannels.
444  using ChanToWireMap = std::unordered_map<raw::ChannelID_t,const recob::Wire*>;
445 
446  ChanToWireMap channelToWireMap;
447 
448  // We will use the presence of a RawDigit as an indicator of a good channel... So
449  // we want a mapping between channel and RawDigit
450  using ChanToRawDigitMap = std::unordered_map<raw::ChannelID_t,const raw::RawDigit*>;
451 
452  ChanToRawDigitMap chanToRawDigitMap;
453 
454  // Look up the list of bad channels
455  art::Handle< std::vector<int>> badChannelHandle;
456  event.getByLabel(fBadChannelProducerLabel, badChannelHandle);
457 
458  // Now start a loop over the individual TPCs to build out the structures for RawDigits and Wires
459  for(size_t tpcID = 0; tpcID < fRawDigitProducerLabelVec.size(); tpcID++)
460  {
461  art::Handle< std::vector<raw::RawDigit> > rawDigitHandle;
462  event.getByLabel(fRawDigitProducerLabelVec[tpcID], rawDigitHandle);
463 
464  art::Handle< std::vector<recob::Wire> > wireHandle;
465  event.getByLabel(fWireProducerLabelVec[tpcID], wireHandle);
466 
467  if (!rawDigitHandle.isValid() || !wireHandle.isValid()) return;
468 
469  for(const auto& wire : *wireHandle) channelToWireMap[wire.Channel()] = &wire;
470 
471  for(const auto& rawDigit : *rawDigitHandle) chanToRawDigitMap[rawDigit.Channel()] = &rawDigit;
472  }
473 
474  // Now we create a data structure to relate hits to their channel ID
475  using ChanToHitVecMap = std::unordered_map<raw::ChannelID_t,std::vector<const recob::Hit*>>;
476 
477  ChanToHitVecMap channelToHitVec;
478 
479  // And now fill it
480  for(const auto& hitLabel : fHitProducerLabelVec)
481  {
482  art::Handle< std::vector<recob::Hit> > hitHandle;
483  event.getByLabel(hitLabel, hitHandle);
484 
485  for(const auto& hit : *hitHandle) channelToHitVec[hit.Channel()].push_back(&hit);
486  }
487 
488  // It is useful to create a mapping between trackID and MCParticle
489  using TrackIDToMCParticleMap = std::unordered_map<int, const simb::MCParticle*>;
490 
491  TrackIDToMCParticleMap trackIDToMCParticleMap;
492 
493  for(const auto& mcParticle : *mcParticleHandle)
494  trackIDToMCParticleMap[mcParticle.TrackId()] = &mcParticle;
495 
496  const lariov::ChannelStatusProvider& chanFilt = art::ServiceHandle<lariov::ChannelStatusService>()->GetProvider();
497 
498  std::vector<int> nSimChannelHitVec = {0,0,0};
499  std::vector<int> nRecobHitVec = {0,0,0};
500  std::vector<int> nFakeHitVec = {0,0,0};
501  std::vector<int> nSimulatedWiresVec = {0,0,0};
502 
503  unsigned int lastwire=-1;
504 
505  auto const clockData = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(event);
506 
507  for(const auto& partToChanInfo : partToChanToTDCToIDEMap)
508  {
509  TrackIDToMCParticleMap::const_iterator trackIDToMCPartItr = trackIDToMCParticleMap.find(partToChanInfo.first);
510 
511  if (trackIDToMCPartItr == trackIDToMCParticleMap.end()) continue;
512 
513  int trackPDGCode = trackIDToMCPartItr->second->PdgCode();
514  std::string processName = trackIDToMCPartItr->second->Process();
515 
516  // Looking for primary muons (e.g. CR Tracks)
517  if (fabs(trackPDGCode) != 13 || processName != "primary") continue;
518 
519  // Recover particle position and angle information
520  Eigen::Vector3f partStartPos(trackIDToMCPartItr->second->Vx(),trackIDToMCPartItr->second->Vy(),trackIDToMCPartItr->second->Vz());
521  Eigen::Vector3f partStartDir(trackIDToMCPartItr->second->Px(),trackIDToMCPartItr->second->Py(),trackIDToMCPartItr->second->Pz());
522 
523  partStartDir.normalize();
524 
525  Eigen::Vector2f partStartDirVecXZ(partStartDir[0],partStartDir[2]);
526 
527  partStartDirVecXZ.normalize();
528 
529  // Assuming the SimChannels contain position information (currently not true for WC produced SimChannels)
530  // then we want to keep a running position
531  std::vector<Eigen::Vector3f> lastPositionVec = {partStartPos,partStartPos,partStartPos};
532 
533  for(const auto& chanToTDCToIDEMap : partToChanInfo.second)
534  {
535  // skip bad channels
536  if (fUseBadChannelDB)
537  {
538  // This is the "correct" way to check and remove bad channels...
539  if( chanFilt.Status(chanToTDCToIDEMap.first) < fMinAllowedChanStatus)
540  {
541  std::vector<geo::WireID> wids = fGeometry->ChannelToWire(chanToTDCToIDEMap.first);
542  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;
543  continue;
544  }
545  }
546 
547  // Was a list made available?
548  // If so then we try that
549  if (badChannelHandle.isValid())
550  {
551  // Here we query the input list from the wirecell processing
552  std::vector<int>::const_iterator badItr = std::find(badChannelHandle->begin(),badChannelHandle->end(),chanToTDCToIDEMap.first);
553 
554  if (badItr != badChannelHandle->end()) continue;
555  // {
556  // ChanToRawDigitMap::const_iterator rawDigitItr = chanToRawDigitMap.find(chanToTDCToIDEMap.first);
557  //
558  // if (rawDigitItr != chanToRawDigitMap.end())
559  // {
560  // float nSig(3.);
561  // float mean(0.);
562  // float rmsFull(0.);
563  // float rmsTrunc(0.);
564  // int nTrunc(0);
565  //
566  // getTruncatedMeanRMS(rawDigitItr->second->ADCs(), nSig, mean, rmsFull, rmsTrunc, nTrunc);
567  //
568  // std::cout << "--> Rejecting channel: " << chanToTDCToIDEMap.first << " from bad channel list, rms: " << rmsFull << std::endl;
569  // }
570  //
571  // continue;
572  // }
573  }
574 
575  TickTDCIDEVec tdcToIDEVec = chanToTDCToIDEMap.second;
576  float totalElectrons(0.);
577  float maxElectrons(0.);
578  unsigned short maxElectronsTDC(0);
579  int nMatchedWires(0);
580  int nMatchedHits(0);
581 
582  // Make sure the vector is time ordered
583  std::sort(tdcToIDEVec.begin(),tdcToIDEVec.end(),[](const auto& left, const auto& right){return left.first < right.first;});
584 
585  // The below try-catch block may no longer be necessary
586  // Decode the channel and make sure we have a valid one
587  std::vector<geo::WireID> wids = fGeometry->ChannelToWire(chanToTDCToIDEMap.first);
588 
589  // Recover plane and wire in the plane
590  unsigned int plane = wids[0].Plane;
591  unsigned int wire = wids[0].Wire;
592 
593  Eigen::Vector3f avePosition(0.,0.,0.);
594 
595  if(wire!=lastwire) nSimulatedWiresVec[plane]++;
596  lastwire=wire;
597 
598  for(const auto& tdcIdePair : tdcToIDEVec)
599  {
600  const sim::IDE* ide = tdcIdePair.second;
601 
602  if (trackPDGCode != ide->trackID) continue;
603 
604  totalElectrons += ide->numElectrons;
605 
606  if (maxElectrons < ide->numElectrons)
607  {
608  maxElectrons = ide->numElectrons;
609  maxElectronsTDC = tdcIdePair.first;
610  }
611 
612  avePosition += Eigen::Vector3f(ide->x,ide->y,ide->z);
613  }
614 
615  // Get local track direction by using the average position of deposited charge as the current position
616  // and then subtracting the last position
617  avePosition /= float(tdcToIDEVec.size());
618 
619  Eigen::Vector3f partDirVec = avePosition - lastPositionVec[plane];
620 
621  partDirVec.normalize();
622 
623  lastPositionVec[plane] = avePosition;
624 
625  // Now what we want is the projected angle in the XZ plane
626  Eigen::Vector2f projPairDirVec(partDirVec[0],partDirVec[2]);
627 
628  projPairDirVec.normalize();
629 
630  float cosThetaXZ = projPairDirVec[0];
631 
632  totalElectrons = std::min(totalElectrons, float(99900.));
633 
634  fTotalElectronsHistVec[plane]->Fill(totalElectrons, 1.);
635  fMaxElectronsHistVec[plane]->Fill(maxElectrons, 1.);
636 
637  nSimChannelHitVec[plane]++;
638 
639  unsigned short startTDC = tdcToIDEVec.begin()->first;
640  unsigned short stopTDC = tdcToIDEVec.rbegin()->first;
641 
642  // Convert to ticks to get in same units as hits
643  unsigned short startTick = clockData.TPCTDC2Tick(startTDC) + fOffsetVec[plane];
644  unsigned short stopTick = clockData.TPCTDC2Tick(stopTDC) + fOffsetVec[plane];
645  unsigned short maxETick = clockData.TPCTDC2Tick(maxElectronsTDC) + fOffsetVec[plane];
646 
647  fSimNumTDCVec[plane]->Fill(stopTick - startTick, 1.);
648 
649  // Set up to extract the "best" parameters in the event of more than one hit for this pulse train
650  float nElectronsTotalBest(0.);
651  float hitSummedADCBest(0.);
652  float hitIntegralBest(0.);
653  float hitPeakTimeBest(0.);
654  float hitPeakAmpBest(-100.);
655  float hitRMSBest(0.);
656  int hitMultiplicityBest(0);
657  int hitLocalIndexBest(0);
658  float hitGoodnessBest(0.);
659  int hitNumDegreesBest(0);
660  float hitBaselineBest(0.);
661  float hitSnippetLenBest(0.);
662  unsigned short hitStopTickBest(0);
663  unsigned short hitStartTickBest(0);
664 
665  // Start by recovering the Wire associated to this channel
666  ChanToWireMap::const_iterator wireItr = channelToWireMap.find(chanToTDCToIDEMap.first);
667 
668  if (wireItr != channelToWireMap.end())
669  {
670  const recob::Wire::RegionsOfInterest_t& signalROI = wireItr->second->SignalROI();
671  const lar::sparse_vector<float>::datarange_t* wireRangePtr(NULL);
672 
673  // Here we need to match the range of the ROI's on the given Wire with the tick range from the SimChannel
674  for(const auto& range : signalROI.get_ranges())
675  {
676  // #################################################
677  // ### Getting a vector of signals for this wire ###
678  // #################################################
679  //std::vector<float> signal(wire->Signal());
680 
681  const std::vector<float>& signal = range.data();
682  raw::TDCtick_t roiFirstBinTick = range.begin_index();
683  raw::TDCtick_t roiLastBinTick = roiFirstBinTick + signal.size();
684 
685  // If no overlap then go to next
686  if (roiFirstBinTick > stopTick || roiLastBinTick < startTick) continue;
687 
688  wireRangePtr = &range;
689  break;
690  }
691 
692  // Check that we have found the wire range
693  // Note that if we have not matched an ROI then we can't have a hit either so skip search for that...
694  if (wireRangePtr)
695  {
696  const recob::Hit* rejectedHit = 0;
697  const recob::Hit* bestHit = 0;
698 
699  nMatchedWires++;
700 
701  // The next mission is to recover the hits associated to this Wire
702  // The easiest way to do this is to simply look up all the hits on this channel and then match
703  ChanToHitVecMap::iterator hitIter = channelToHitVec.find(chanToTDCToIDEMap.first);
704 
705  if (hitIter != channelToHitVec.end())
706  {
707  // Loop through the hits for this channel and look for matches
708  // In the event of more than one hit associated to the sim channel range, keep only
709  // the best match (assuming the nearby hits are "extra")
710  // Note that assumption breaks down for long pulse trains but worry about that later
711  for(const auto& hit : hitIter->second)
712  {
713  unsigned short hitStartTick = hit->PeakTime() - fSigmaVec[plane] * hit->RMS();
714  unsigned short hitStopTick = hit->PeakTime() + fSigmaVec[plane] * hit->RMS();
715 
716  // If hit is out of range then skip, it is not related to this particle
717  if (hitStartTick > stopTick || hitStopTick < startTick)
718  {
719  nFakeHitVec[plane]++;
720  rejectedHit = hit;
721  continue;
722  }
723 
724  float hitHeight = hit->PeakAmplitude();
725 
726  // Use the hit with the largest pulse height as the "best"
727  if (hitHeight < hitPeakAmpBest) continue;
728 
729  hitPeakAmpBest = hitHeight;
730  bestHit = hit;
731  hitStartTickBest = hitStartTick;
732  hitStopTickBest = hitStopTick;
733  }
734 
735  // Find a match?
736  if (bestHit)
737  {
738  nElectronsTotalBest = 0.;
739  hitPeakTimeBest = bestHit->PeakTime();
740  hitIntegralBest = bestHit->Integral();
741  hitSummedADCBest = bestHit->SummedADC();
742  hitRMSBest = bestHit->RMS();
743  hitMultiplicityBest = bestHit->Multiplicity();
744  hitLocalIndexBest = bestHit->LocalIndex();
745  hitGoodnessBest = bestHit->GoodnessOfFit();
746  hitNumDegreesBest = bestHit->DegreesOfFreedom();
747  hitSnippetLenBest = bestHit->EndTick() - bestHit->StartTick();
748  hitBaselineBest = 0.; // To do...
749 
750  nMatchedHits++;
751 
752  // Get the number of electrons
753  for(int tickIdx = 0; tickIdx < int(tdcToIDEVec.size()); tickIdx++)
754  {
755  // We might be done?
756  if (tickIdx > hitStopTickBest - hitStartTickBest) break;
757 
758  // Convert to TDC
759  unsigned short hitTDC = clockData.TPCTick2TDC(hitStartTickBest + tickIdx - fOffsetVec[plane]);
760 
761  if (hitTDC >= tdcToIDEVec[tickIdx].first)
762  {
763  if (tdcToIDEVec[tickIdx].second->trackID == trackPDGCode) nElectronsTotalBest += tdcToIDEVec[tickIdx].second->numElectrons;
764  }
765  }
766  }
767 
768  if (nMatchedHits > 0)
769  {
770  float chgRatioADC = hitSummedADCBest > 0. ? totalElectrons / hitSummedADCBest : 0.;
771  float chgRatioInt = hitIntegralBest > 0. ? totalElectrons / hitIntegralBest : 0.;
772 
773  fHitSumADCVec[plane]->Fill(hitSummedADCBest, 1.);
774  fHitIntegralHistVec[plane]->Fill(hitIntegralBest, 1.);
775  fSimDivHitChgVec[plane]->Fill(chgRatioADC, 1.);
776  fSimDivHitChg1Vec[plane]->Fill(chgRatioInt, 1.);
777  fHitVsSimChgVec[plane]->Fill(std::min(hitSummedADCBest,float(999.)), totalElectrons, 1.);
778  fHitVsSimIntVec[plane]->Fill(std::min(hitIntegralBest,float(999.)), totalElectrons, 1.);
779  fToteVHitEIntVec[plane]->Fill(std::min(totalElectrons,float(99999.)),std::min(nElectronsTotalBest,float(99999.)),1.);
780  fHitPulseHeightVec[plane]->Fill(std::min(hitPeakAmpBest,float(149.5)), 1.);
781  fHitPulseWidthVec[plane]->Fill(std::min(hitRMSBest,float(19.8)), 1.);
782  fHitElectronsVec[plane]->Fill(nElectronsTotalBest, 1.);
783  fHitNumTDCVec[plane]->Fill(std::min(float(hitStopTickBest - hitStartTickBest),float(99.5)), 1.);
784  fSnippetLenVec[plane]->Fill(std::min(hitSnippetLenBest, float(99.5)), 1.);
785  fDeltaMidTDCVec[plane]->Fill(hitPeakTimeBest - maxETick, 1.);
786 
787  nRecobHitVec[plane]++;
788 
789  }
790  else if (rejectedHit)
791  {
792  unsigned short hitStartTick = rejectedHit->PeakTime() - fSigmaVec[plane] * rejectedHit->RMS();
793  unsigned short hitStopTick = rejectedHit->PeakTime() + fSigmaVec[plane] * rejectedHit->RMS();
794 
795  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;
796  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;
797  }
798  else
799  {
800  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;
801  }
802  }
803  }
804  }
805 
806  fWireEfficVec.at(plane)->Fill(totalElectrons, std::min(nMatchedWires,1), 1.);
807  fWireEfficPHVec.at(plane)->Fill(maxElectrons, std::min(nMatchedWires,1), 1.);
808 
809  float matchHit = std::min(nMatchedHits,1);
810  float snippetLen = std::min(float(stopTick - startTick),float(99.5));
811 
812  fNMatchedHitVec[plane]->Fill(nMatchedHits, 1.);
813  fHitEfficVec[plane]->Fill(totalElectrons, matchHit, 1.);
814  fHitEfficPHVec[plane]->Fill(maxElectrons, matchHit, 1.);
815  fHitEfficXZVec[plane]->Fill(cosThetaXZ, matchHit, 1.);
816  fCosXZvRMSVec[plane]->Fill(cosThetaXZ, snippetLen, 1.);
817  fHitEfficRMSVec[plane]->Fill(snippetLen, matchHit, 1.);
818 
819  fHitENEvXZVec[plane]->Fill(cosThetaXZ, totalElectrons, matchHit);
820 
821  // Store tuple variables
822  fTPCVec.push_back(wids[0].TPC);
823  fCryoVec.push_back(wids[0].Cryostat);
824  fPlaneVec.push_back(wids[0].Plane);
825  fWireVec.push_back(wids[0].Wire);
826 
827  fTotalElectronsVec.push_back(totalElectrons);
828  fMaxElectronsVec.push_back(maxElectrons);
829  fStartTickVec.push_back(startTick);
830  fStopTickVec.push_back(stopTick);
831  fMaxETickVec.push_back(maxETick);
832  fPartDirX.push_back(partDirVec[0]);
833  fPartDirY.push_back(partDirVec[1]);
834  fPartDirZ.push_back(partDirVec[2]);
835 
836  fNMatchedWires.push_back(nMatchedWires);
837  fNMatchedHits.push_back(nMatchedHits);
838 
839  fHitPeakTimeVec.push_back(hitPeakTimeBest);
840  fHitPeakAmpVec.push_back(hitPeakAmpBest);
841  fHitPeakRMSVec.push_back(hitRMSBest);
842  fHitBaselinevec.push_back(hitBaselineBest);
843  fHitSummedADCVec.push_back(hitSummedADCBest);
844  fHitIntegralVec.push_back(hitIntegralBest);
845  fHitStartTickVec.push_back(hitStartTickBest);
846  fHitStopTickVec.push_back(hitStopTickBest);
847  fHitMultiplicityVec.push_back(hitMultiplicityBest);
848  fHitLocalIndexVec.push_back(hitLocalIndexBest);
849  fHitGoodnessVec.push_back(hitGoodnessBest);
850  fNumDegreesVec.push_back(hitNumDegreesBest);
851  }
852  }
853 
854  for(size_t idx = 0; idx < fGeometry->Nplanes();idx++)
855  {
856  if (nSimChannelHitVec[idx] > 10)
857  {
858  float hitEfficiency = float(nRecobHitVec[idx]) / float(nSimChannelHitVec[idx]);
859 
860  fNSimChannelHitsVec[idx]->Fill(std::min(nSimChannelHitVec[idx],1999),1.);
861  fNRecobHitVec[idx]->Fill(std::min(nRecobHitVec[idx],1999), 1.);
862  fNFakeHitVec[idx]->Fill(nFakeHitVec[idx]/(float)nSimulatedWiresVec[idx],1.);
863  fHitEfficiencyVec[idx]->Fill(hitEfficiency, 1.);
864  }
865  }
866 
867  return;
868 }
869 
870 // Useful for normalizing histograms
871 void TrackHitEfficiencyAnalysis::endJob(int numEvents)
872 {
873  return;
874 }
875 
876 DEFINE_ART_CLASS_TOOL(TrackHitEfficiencyAnalysis)
877 }
short int LocalIndex() const
How well do we believe we know this hit?
Definition: Hit.h:227
TrackID_t trackID
Geant4 supplied track ID.
Definition: SimChannel.h:116
float z
z position of ionization [cm]
Definition: SimChannel.h:121
Utilities related to art service access.
BEGIN_PROLOG true icarus_rawdigitfilter FilterTools FilterPlane1 Plane
geo::WireID WireID() const
Definition: Hit.h:233
float RMS() const
RMS of the hit shape, in tick units.
Definition: Hit.h:220
Declaration of signal hit object.
walls no right
Definition: selectors.fcl:105
int DegreesOfFreedom() const
Definition: Hit.h:229
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.
const range_list_t & get_ranges() const
Returns the internal list of non-void ranges.
float Integral() const
Integral under the calibrated signal waveform of the hit, in tick x ADC units.
Definition: Hit.h:224
WireID_t Wire
Index of the wire within its plane.
Definition: geo_types.h:580
process_name hit
Definition: cheaterreco.fcl:51
float GoodnessOfFit() const
Degrees of freedom in the determination of the hit signal shape (-1 by default)
Definition: Hit.h:228
short int Multiplicity() const
How many hits could this one be shared with.
Definition: Hit.h:226
float x
x position of ionization [cm]
Definition: SimChannel.h:119
BEGIN_PROLOG TPC
float PeakAmplitude() const
The estimated amplitude of the hit at its peak, in ADC units.
Definition: Hit.h:221
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.
Definition: RawTypes.h:25
void endJob(int numEvents) override
Interface for method to executve at the end of run processing.
unsigned int Nplanes(unsigned int tpc=0, unsigned int cstat=0) const
Returns the total number of wire planes in the specified TPC.
void fillHistograms(const art::Event &) const override
Interface for filling histograms.
Ionization at a point of the TPC sensitive volume.
Definition: SimChannel.h:86
void initializeTuple(TTree *) override
Interface for initializing the tuple variables.
virtual Status_t Status(raw::ChannelID_t channel) const
Returns a status integer with arbitrary meaning.
walls no left
Definition: selectors.fcl:105
Class providing information about the quality of channels.
PlaneID_t Plane
Index of the plane within its TPC.
Definition: geo_types.h:493
Description of geometry of one entire detector.
raw::TDCtick_t StartTick() const
Initial tdc tick for hit.
Definition: Hit.h:216
tuple dir
Definition: dropbox.py:28
float y
y position of ionization [cm]
Definition: SimChannel.h:120
raw::TDCtick_t EndTick() const
Final tdc tick for hit.
Definition: Hit.h:217
float PeakTime() const
Time of the signal peak, in tick units.
Definition: Hit.h:218
std::string to_string(WindowPattern const &pattern)
object containing MC truth information necessary for making RawDigits and doing back tracking ...
Interface for experiment-specific channel quality info provider.
Range class, with range and data.
float SummedADC() const
The sum of calibrated ADC counts of the hit (0. by default)
Definition: Hit.h:223
Declaration of basic channel signal object.
2D representation of charge deposited in the TDC/wire plane
Definition: Hit.h:48
art::ServiceHandle< art::TFileService > tfs
unsigned int ChannelID_t
Type representing the ID of a readout channel.
Definition: RawTypes.h:28
TPCID_t TPC
Index of the TPC within its cryostat.
Definition: geo_types.h:406
Interface for experiment-specific service for channel quality info.
art framework interface to geometry description
BEGIN_PROLOG could also be cout
float numElectrons
number of electrons at the readout for this track ID and time
Definition: SimChannel.h:117