All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
sbndcode/sbndcode/Calibration/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"
11 
13 #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 
35 namespace TrackHitEfficiencyAnalysis
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;
113  art::InputTag fMCParticleProducerLabel;
114  art::InputTag fSimChannelProducerLabel;
115  art::InputTag fBadChannelProducerLabel;
116  bool fUseBadChannelDB;
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
121 
122  // Pointers to the histograms we'll create.
123  std::vector<TH1F*> fTotalElectronsHistVec;
124  std::vector<TH1F*> fMaxElectronsHistVec;
125  std::vector<TH1F*> fHitElectronsVec;
126  std::vector<TH1F*> fHitSumADCVec;
127  std::vector<TH1F*> fHitIntegralHistVec;
128  std::vector<TH1F*> fHitPulseHeightVec;
129  std::vector<TH1F*> fHitPulseWidthVec;
130  std::vector<TH1F*> fSimNumTDCVec;
131  std::vector<TH1F*> fHitNumTDCVec;
132  std::vector<TH1F*> fSnippetLenVec;
133  std::vector<TH1F*> fNMatchedHitVec;
134  std::vector<TH1F*> fDeltaMidTDCVec;
135  std::vector<TProfile*> fWireEfficVec;
136  std::vector<TProfile*> fWireEfficPHVec;
137  std::vector<TProfile*> fHitEfficVec;
138  std::vector<TProfile*> fHitEfficPHVec;
139  std::vector<TProfile*> fHitEfficXZVec;
140  std::vector<TProfile*> fHitEfficRMSVec;
141  std::vector<TProfile*> fCosXZvRMSVec;
142  std::vector<TProfile2D*> fHitENEvXZVec;
143  std::vector<TH1F*> fSimDivHitChgVec;
144  std::vector<TH1F*> fSimDivHitChg1Vec;
145  std::vector<TH2F*> fHitVsSimChgVec;
146  std::vector<TH2F*> fHitVsSimIntVec;
147  std::vector<TH2F*> fToteVHitEIntVec;
148 
149  std::vector<TH1F*> fNSimChannelHitsVec;
150  std::vector<TH1F*> fNRecobHitVec;
151  //std::vector<TH1F*> fNRejectedHitVec;
152  std::vector<TH1F*> fHitEfficiencyVec;
153 
154  std::vector<TH1F*> fNFakeHitVec;
155 
156  // TTree variables
157  mutable TTree* fTree;
158 
159  mutable std::vector<int> fTPCVec;
160  mutable std::vector<int> fCryoVec;
161  mutable std::vector<int> fPlaneVec;
162  mutable std::vector<int> fWireVec;
163 
164  mutable std::vector<float> fTotalElectronsVec;
165  mutable std::vector<float> fMaxElectronsVec;
166  mutable std::vector<int> fStartTickVec;
167  mutable std::vector<int> fStopTickVec;
168  mutable std::vector<int> fMaxETickVec;
169  mutable std::vector<float> fPartDirX;
170  mutable std::vector<float> fPartDirY;
171  mutable std::vector<float> fPartDirZ;
172 
173  mutable std::vector<int> fNMatchedWires;
174  mutable std::vector<int> fNMatchedHits;
175 
176  mutable std::vector<float> fHitPeakTimeVec;
177  mutable std::vector<float> fHitPeakAmpVec;
178  mutable std::vector<float> fHitPeakRMSVec;
179  mutable std::vector<float> fHitBaselinevec;
180  mutable std::vector<float> fHitSummedADCVec;
181  mutable std::vector<float> fHitIntegralVec;
182  mutable std::vector<int> fHitStartTickVec;
183  mutable std::vector<int> fHitStopTickVec;
184  mutable std::vector<int> fHitMultiplicityVec;
185  mutable std::vector<int> fHitLocalIndexVec;
186  mutable std::vector<float> fHitGoodnessVec;
187  mutable std::vector<int> fNumDegreesVec;
188 
189  // Useful services, keep copies for now (we can update during begin run periods)
190  const geo::GeometryCore* fGeometry; ///< pointer to Geometry service
191 };
192 
193 //----------------------------------------------------------------------------
194 /// Constructor.
195 ///
196 /// Arguments:
197 ///
198 /// pset - Fcl parameters.
199 ///
200 TrackHitEfficiencyAnalysis::TrackHitEfficiencyAnalysis(fhicl::ParameterSet const & pset) : fTree(nullptr)
201 {
202  fGeometry = lar::providerFrom<geo::Geometry>();
203 
204  configure(pset);
205 
206  // Report.
207  mf::LogInfo("TrackHitEfficiencyAnalysis") << "TrackHitEfficiencyAnalysis configured\n";
208 }
209 
210 //----------------------------------------------------------------------------
211 /// Destructor.
212 TrackHitEfficiencyAnalysis::~TrackHitEfficiencyAnalysis()
213 {}
214 
215 //----------------------------------------------------------------------------
216 /// Reconfigure method.
217 ///
218 /// Arguments:
219 ///
220 /// pset - Fcl parameter set.
221 ///
222 void TrackHitEfficiencyAnalysis::configure(fhicl::ParameterSet const & pset)
223 {
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"});
227  fMCParticleProducerLabel = pset.get< art::InputTag >("MCParticleLabel", "largeant");
228  fSimChannelProducerLabel = pset.get< art::InputTag >("SimChannelLabel", "largeant");
229  fBadChannelProducerLabel = pset.get< art::InputTag >("BadChannelLabel", "simnfspl1:badchannels");
230  fUseBadChannelDB = pset.get< bool >("UseBadChannelDB", true);
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.});
234  fMinAllowedChanStatus = pset.get< int >("MinAllowedChannelStatus");
235 }
236 
237 //----------------------------------------------------------------------------
238 /// Begin job method.
239 void TrackHitEfficiencyAnalysis::initializeHists(art::ServiceHandle<art::TFileService>& tfs, const std::string& dirName)
240 {
241  // Make a directory for these histograms
242  art::TFileDirectory dir = tfs->mkdir(dirName.c_str());
243 
247  fHitSumADCVec.resize(fGeometry->Nplanes());
251  fSimNumTDCVec.resize(fGeometry->Nplanes());
252  fHitNumTDCVec.resize(fGeometry->Nplanes());
253  fSnippetLenVec.resize(fGeometry->Nplanes());
254  fNMatchedHitVec.resize(fGeometry->Nplanes());
255  fDeltaMidTDCVec.resize(fGeometry->Nplanes());
256  fHitVsSimChgVec.resize(fGeometry->Nplanes());
257  fHitVsSimIntVec.resize(fGeometry->Nplanes());
258  fCosXZvRMSVec.resize(fGeometry->Nplanes());
260  fHitENEvXZVec.resize(fGeometry->Nplanes());
262  fNRecobHitVec.resize(fGeometry->Nplanes());
263  fNFakeHitVec.resize(fGeometry->Nplanes());
267 
268  fWireEfficVec.resize(fGeometry->Nplanes());
269  fWireEfficPHVec.resize(fGeometry->Nplanes());
270 
271  fHitEfficVec.resize(fGeometry->Nplanes());
272  fHitEfficPHVec.resize(fGeometry->Nplanes());
273  fHitEfficXZVec.resize(fGeometry->Nplanes());
274  fHitEfficRMSVec.resize(fGeometry->Nplanes());
275 
276  for(size_t plane = 0; plane < fGeometry->Nplanes(); plane++)
277  {
278  fTotalElectronsHistVec.at(plane) = dir.make<TH1F>(("TotalElecs" + std::to_string(plane)).c_str(), ";Total # electrons", 250, 0., 100000.);
279  fMaxElectronsHistVec.at(plane) = dir.make<TH1F>(("MaxElecs" + std::to_string(plane)).c_str(), ";Max # electrons", 250, 0., 20000.);
280  fHitElectronsVec.at(plane) = dir.make<TH1F>(("HitElecs" + std::to_string(plane)).c_str(), ";# e- in Hit Range", 250, 0., 100000.);
281  fHitSumADCVec.at(plane) = dir.make<TH1F>(("SumADC" + std::to_string(plane)).c_str(), "Hit Sum ADC", 200, 0., 1000.);
282  fHitIntegralHistVec.at(plane) = dir.make<TH1F>(("Integral" + std::to_string(plane)).c_str(), "Hit Integral", 200, 0., 1000.);
283  fHitPulseHeightVec.at(plane) = dir.make<TH1F>(("PulseHeight" + std::to_string(plane)).c_str(), "Hit PH (ADC)", 200, 0., 150.);
284  fHitPulseWidthVec.at(plane) = dir.make<TH1F>(("PulseWidth" + std::to_string(plane)).c_str(), "Hit PulseWidth;Hit RMS", 100, 0., 20.);
285  fSimNumTDCVec.at(plane) = dir.make<TH1F>(("SimNumTDC" + std::to_string(plane)).c_str(), ";TDC ticks", 100, 0., 100.);
286  fHitNumTDCVec.at(plane) = dir.make<TH1F>(("HitNumTDC" + std::to_string(plane)).c_str(), ";ticks", 100, 0., 100.);
287  fSnippetLenVec.at(plane) = dir.make<TH1F>(("SnippetLen" + std::to_string(plane)).c_str(), ";ticks", 100, 0., 100.);
288  fNMatchedHitVec.at(plane) = dir.make<TH1F>(("NMatched" + std::to_string(plane)).c_str(), ";# hits", 20, 0., 20.);
289  fDeltaMidTDCVec.at(plane) = dir.make<TH1F>(("DeltaMid" + std::to_string(plane)).c_str(), ";# hits", 50, -25., 25.);
290  fNSimChannelHitsVec.at(plane) = dir.make<TH1F>(("NSimChan" + std::to_string(plane)).c_str(), ";# hits", 100, 0., 1200.);
291  fNRecobHitVec.at(plane) = dir.make<TH1F>(("NRecobHit" + std::to_string(plane)).c_str(), ";# hits", 100, 0., 1200.);
292  fNFakeHitVec.at(plane) = dir.make<TH1F>(("NFakeHit" + std::to_string(plane)).c_str(), ";# hits", 100, 0., 50.);
293  fHitEfficiencyVec.at(plane) = dir.make<TH1F>(("PlnEffic" + std::to_string(plane)).c_str(), ";# hits", 101, 0., 1.01);
294  fSimDivHitChgVec.at(plane) = dir.make<TH1F>(("SimDivHit" + std::to_string(plane)).c_str(), ";# e / SummedADC", 200, 0., 200.);
295  fSimDivHitChg1Vec.at(plane) = dir.make<TH1F>(("SimDivHit1" + std::to_string(plane)).c_str(), ";# e / Integral", 200, 0., 200.);
296 
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.);
300 
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.);
303 
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.);
311  }
312 
313  return;
314 }
315 
316 void TrackHitEfficiencyAnalysis::initializeTuple(TTree* tree)
317 {
318  fTree = tree;
319 
320  fTree->Branch("CryostataVec", "std::vector<int>", &fCryoVec);
321  fTree->Branch("TPCVec", "std::vector<int>", &fTPCVec);
322  fTree->Branch("PlaneVec", "std::vector<int>", &fPlaneVec);
323  fTree->Branch("WireVec", "std::vector<int>", &fWireVec);
324 
325  fTree->Branch("TotalElectronsVec", "std::vector<float>", &fTotalElectronsVec);
326  fTree->Branch("MaxElectronsVec", "std::vector<float>", &fMaxElectronsVec);
327  fTree->Branch("StartTick", "std::vector<int>", &fStartTickVec);
328  fTree->Branch("StopTick", "std::vector<int>", &fStopTickVec);
329  fTree->Branch("MaxETick", "std::vector<int>", &fMaxETickVec);
330  fTree->Branch("PartDirX", "std::vector<float>", &fPartDirX);
331  fTree->Branch("PartDirY", "std::vector<float>", &fPartDirY);
332  fTree->Branch("PartDirZ", "std::vector<float>", &fPartDirZ);
333 
334  fTree->Branch("NMatchedWires", "std::vector<int>", &fNMatchedWires);
335  fTree->Branch("NMatchedHits", "std::vector<int>", &fNMatchedHits);
336 
337  fTree->Branch("HitPeakTimeVec", "std::vector<float>", &fHitPeakTimeVec);
338  fTree->Branch("HitPeakAmpVec", "std::vector<float>", &fHitPeakAmpVec);
339  fTree->Branch("HitPeakRMSVec", "std::vector<float>", &fHitPeakRMSVec);
340  fTree->Branch("HitBaselineVec", "std::vector<float>", &fHitBaselinevec);
341  fTree->Branch("HitSummedADCVec", "std::vector<float>", &fHitSummedADCVec);
342  fTree->Branch("HitIntegralVec", "std::vector<float>", &fHitIntegralVec);
343  fTree->Branch("HitStartTickVec", "std::vector<int>", &fHitStartTickVec);
344  fTree->Branch("HitStopTickVec", "std::vector<int>", &fHitStopTickVec);
345  fTree->Branch("HitMultiplicity", "std::vector<int>", &fHitMultiplicityVec);
346  fTree->Branch("HitLocalIndex", "std::vector<int>", &fHitLocalIndexVec);
347  fTree->Branch("HitGoodness", "std::vector<float>", &fHitGoodnessVec);
348  fTree->Branch("HitNumDegrees", "std::vector<int>", &fNumDegreesVec);
349 
350  clear();
351 
352  return;
353 }
354 
355 void TrackHitEfficiencyAnalysis::clear() const
356 {
357  fTPCVec.clear();
358  fCryoVec.clear();
359  fPlaneVec.clear();
360  fWireVec.clear();
361 
362  fTotalElectronsVec.clear();
363  fMaxElectronsVec.clear();
364  fStartTickVec.clear();
365  fStopTickVec.clear();
366  fMaxETickVec.clear();
367  fPartDirX.clear();
368  fPartDirY.clear();
369  fPartDirZ.clear();
370 
371  fNMatchedWires.clear();
372  fNMatchedHits.clear();
373 
374  fHitPeakTimeVec.clear();
375  fHitPeakAmpVec.clear();
376  fHitPeakRMSVec.clear();
377  fHitBaselinevec.clear();
378  fHitSummedADCVec.clear();
379  fHitIntegralVec.clear();
380  fHitStartTickVec.clear();
381  fHitStopTickVec.clear();
382  fHitMultiplicityVec.clear();
383  fHitLocalIndexVec.clear();
384  fHitGoodnessVec.clear();
385  fNumDegreesVec.clear();
386 
387  return;
388 }
389 
390 void TrackHitEfficiencyAnalysis::fillHistograms(const art::Event& event) const
391 {
392  std::cout << " filling histos " << std::endl;
393  // Basic assumption is that the producer label vecs for RawDigits and Wire data are
394  // all the same length and in the same order. Here we just check for length
395  if (fRawDigitProducerLabelVec.size() != fWireProducerLabelVec.size()) return;
396  std::cout<<fRawDigitProducerLabelVec.size()<<std::endl;
397  // Always clear the tuple
398  clear();
399 
400  art::Handle< std::vector<sim::SimChannel>> simChannelHandle;
401  event.getByLabel(fSimChannelProducerLabel, simChannelHandle);
402 
403  art::Handle< std::vector<simb::MCParticle>> mcParticleHandle;
404  event.getByLabel(fMCParticleProducerLabel, mcParticleHandle);
405 
406  // If there is no sim channel informaton then exit
407  if (!simChannelHandle.isValid() || simChannelHandle->empty() || !mcParticleHandle.isValid())
408  {
409  std::cout<<"No sim channel information"<<std::endl;
410  return;
411  }
412  std::cout<<"There is sim channel information"<<std::endl;
413 
414 
415  // There are several things going on here... for each channel we have particles (track id's) depositing energy in a range to ticks
416  // So... for each channel we want to build a structure that relates particles to tdc ranges and deposited energy (or electrons)
417  // Here is a complicated structure:
418  using TDCToIDEMap = std::map<unsigned short, sim::IDE>; // We need this one in order
419  using ChanToTDCToIDEMap = std::map<raw::ChannelID_t, TDCToIDEMap>;
420  using PartToChanToTDCToIDEMap = std::unordered_map<int, ChanToTDCToIDEMap>;
421 
422  PartToChanToTDCToIDEMap partToChanToTDCToIDEMap;
423 
424  // Build out the above data structure
425  for(const auto& simChannel : *simChannelHandle)
426  {
427  for(const auto& tdcide : simChannel.TDCIDEMap())
428  {
429  for(const auto& ide : tdcide.second) partToChanToTDCToIDEMap[ide.trackID][simChannel.Channel()][tdcide.first] = ide;
430  }
431  }
432 
433  // what needs to be done?
434  // First we define a straightforward channel to Wire map so we can look up a given
435  // channel's Wire data as we loop over SimChannels.
436  using ChanToWireMap = std::unordered_map<raw::ChannelID_t,const recob::Wire*>;
437 
438  ChanToWireMap channelToWireMap;
439 
440  // We will use the presence of a RawDigit as an indicator of a good channel... So
441  // we want a mapping between channel and RawDigit
442  using ChanToRawDigitMap = std::unordered_map<raw::ChannelID_t,const raw::RawDigit*>;
443 
444  ChanToRawDigitMap chanToRawDigitMap;
445 
446  // Look up the list of bad channels
447  art::Handle< std::vector<int>> badChannelHandle;
448  event.getByLabel(fBadChannelProducerLabel, badChannelHandle);
449 
450  // std::cout<<"Pre TPC loop (L455ish)"<<std::endl;
451 
452  // Now start a loop over the individual TPCs to build out the structures for RawDigits and Wires
453  for(size_t tpcID = 0; tpcID < fRawDigitProducerLabelVec.size(); tpcID++)
454  {
455  //std::cout<<"starting TPC loop"<<std::endl;
456  art::Handle< std::vector<raw::RawDigit> > rawDigitHandle;
457  event.getByLabel(fRawDigitProducerLabelVec[tpcID], rawDigitHandle);
458 
459  //std::cout << "rawDigitHandle: " << rawDigitHandle << std::endl;
460  art::Handle< std::vector<recob::Wire> > wireHandle;;
461  event.getByLabel(fWireProducerLabelVec[tpcID], wireHandle);
462 
463  //std::cout << "wireHandle" << wireHandle << std::endl;
464 
465  if (!rawDigitHandle.isValid() || !wireHandle.isValid())
466  {
467  // std::cout<<"Inside TPC loop if statement (L470ish)"<<std::endl;
468  return;
469  }
470  for(const auto& wire : *wireHandle) channelToWireMap[wire.Channel()] = &wire;
471 
472  for(const auto& rawDigit : *rawDigitHandle) chanToRawDigitMap[rawDigit.Channel()] = &rawDigit;
473  }
474  // std::cout<<"Post TPC loop (L470ish)"<<std::endl;
475 
476  // Now we create a data structure to relate hits to their channel ID
477  using ChanToHitVecMap = std::unordered_map<raw::ChannelID_t,std::vector<const recob::Hit*>>;
478 
479  ChanToHitVecMap channelToHitVec;
480 
481  // And now fill it
482  for(const auto& hitLabel : fHitProducerLabelVec)
483  {
484  art::Handle< std::vector<recob::Hit> > hitHandle;
485  event.getByLabel(hitLabel, hitHandle);
486 
487  for(const auto& hit : *hitHandle) channelToHitVec[hit.Channel()].push_back(&hit);
488  }
489 
490  // It is useful to create a mapping between trackID and MCParticle
491  using TrackIDToMCParticleMap = std::unordered_map<int, const simb::MCParticle*>;
492 
493  TrackIDToMCParticleMap trackIDToMCParticleMap;
494 
495  for(const auto& mcParticle : *mcParticleHandle)
496  trackIDToMCParticleMap[mcParticle.TrackId()] = &mcParticle;
497 
498  const lariov::ChannelStatusProvider& chanFilt = art::ServiceHandle<lariov::ChannelStatusService>()->GetProvider();
499  auto const clockData = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(event);
500 
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};
505 
506  unsigned int lastwire=-1;
507  //std::cout<<"Just before eternal loop of loops (L502)"<<std::endl;
508 
509  for(const auto& partToChanInfo : partToChanToTDCToIDEMap)
510  {
511  //std::cout<<"Starting eternal loop of loops"<<std::endl;
512  TrackIDToMCParticleMap::const_iterator trackIDToMCPartItr = trackIDToMCParticleMap.find(partToChanInfo.first);
513 
514  if (trackIDToMCPartItr == trackIDToMCParticleMap.end()) continue;
515 
516  int trackPDGCode = trackIDToMCPartItr->second->PdgCode();
517  std::string processName = trackIDToMCPartItr->second->Process();
518 
519  // Looking for primary muons (e.g. CR Tracks)
520  if (fabs(trackPDGCode) != 13 || processName != "primary") continue;
521 
522  // Recover particle position and angle information
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());
525 
526  partStartDir.normalize();
527 
528  Eigen::Vector2f partStartDirVecXZ(partStartDir[0],partStartDir[2]);
529 
530  partStartDirVecXZ.normalize();
531 
532  // Assuming the SimChannels contain position information (currently not true for WC produced SimChannels)
533  // then we want to keep a running position
534  std::vector<Eigen::Vector3f> lastPositionVec = {partStartPos,partStartPos,partStartPos};
535 
536  for(const auto& chanToTDCToIDEMap : partToChanInfo.second)
537  {
538  // skip bad channels
539  if (fUseBadChannelDB)
540  {
541  // This is the "correct" way to check and remove bad channels...
542  if( chanFilt.Status(chanToTDCToIDEMap.first) < fMinAllowedChanStatus)
543  {
544  std::vector<geo::WireID> wids = fGeometry->ChannelToWire(chanToTDCToIDEMap.first);
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;
546  continue;
547  }
548  }
549 
550  // Was a list made available?
551  // If so then we try that
552  if (badChannelHandle.isValid())
553  {
554  // Here we query the input list from the wirecell processing
555  std::vector<int>::const_iterator badItr = std::find(badChannelHandle->begin(),badChannelHandle->end(),chanToTDCToIDEMap.first);
556 
557  if (badItr != badChannelHandle->end()) continue;
558  // {
559  // ChanToRawDigitMap::const_iterator rawDigitItr = chanToRawDigitMap.find(chanToTDCToIDEMap.first);
560  //
561  // if (rawDigitItr != chanToRawDigitMap.end())
562  // {
563  // float nSig(3.);
564  // float mean(0.);
565  // float rmsFull(0.);
566  // float rmsTrunc(0.);
567  // int nTrunc(0);
568  //
569  // getTruncatedMeanRMS(rawDigitItr->second->ADCs(), nSig, mean, rmsFull, rmsTrunc, nTrunc);
570  //
571  // std::cout << "--> Rejecting channel: " << chanToTDCToIDEMap.first << " from bad channel list, rms: " << rmsFull << std::endl;
572  // }
573  //
574  // continue;
575  // }
576  }
577 
578  TDCToIDEMap tdcToIDEMap = chanToTDCToIDEMap.second;
579  float totalElectrons(0.);
580  float maxElectrons(0.);
581  unsigned short maxElectronsTDC(0);
582  int nMatchedWires(0);
583  int nMatchedHits(0);
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& ideVal : tdcToIDEMap)
599  {
600  totalElectrons += ideVal.second.numElectrons;
601 
602  if (maxElectrons < ideVal.second.numElectrons)
603  {
604  maxElectrons = ideVal.second.numElectrons;
605  maxElectronsTDC = ideVal.first;
606  }
607 
608  avePosition += Eigen::Vector3f(ideVal.second.x,ideVal.second.y,ideVal.second.z);
609  }
610 
611  // Get local track direction by using the average position of deposited charge as the current position
612  // and then subtracting the last position
613  avePosition /= float(tdcToIDEMap.size());
614 
615  Eigen::Vector3f partDirVec = avePosition - lastPositionVec[plane];
616 
617  partDirVec.normalize();
618 
619  lastPositionVec[plane] = avePosition;
620 
621  // Now what we want is the projected angle in the XZ plane
622  Eigen::Vector2f projPairDirVec(partDirVec[0],partDirVec[2]);
623 
624  projPairDirVec.normalize();
625 
626  float cosThetaXZ = projPairDirVec[0];
627 
628  totalElectrons = std::min(totalElectrons, float(99900.));
629 
630  fTotalElectronsHistVec[plane]->Fill(totalElectrons, 1.);
631  fMaxElectronsHistVec[plane]->Fill(maxElectrons, 1.);
632 
633  nSimChannelHitVec[plane]++;
634 
635  unsigned short startTDC = tdcToIDEMap.begin()->first;
636  unsigned short stopTDC = tdcToIDEMap.rbegin()->first;
637 
638  // Convert to ticks to get in same units as hits
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];
642 
643  fSimNumTDCVec[plane]->Fill(stopTick - startTick, 1.);
644 
645  // Set up to extract the "best" parameters in the event of more than one hit for this pulse train
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);
660 
661  // Start by recovering the Wire associated to this channel
662  ChanToWireMap::const_iterator wireItr = channelToWireMap.find(chanToTDCToIDEMap.first);
663 
664  if (wireItr != channelToWireMap.end())
665  {
666  const recob::Wire::RegionsOfInterest_t& signalROI = wireItr->second->SignalROI();
667  const lar::sparse_vector<float>::datarange_t* wireRangePtr(NULL);
668 
669  // Here we need to match the range of the ROI's on the given Wire with the tick range from the SimChannel
670  for(const auto& range : signalROI.get_ranges())
671  {
672  // #################################################
673  // ### Getting a vector of signals for this wire ###
674  // #################################################
675  //std::vector<float> signal(wire->Signal());
676 
677  const std::vector<float>& signal = range.data();
678  raw::TDCtick_t roiFirstBinTick = range.begin_index();
679  raw::TDCtick_t roiLastBinTick = roiFirstBinTick + signal.size();
680 
681  // If no overlap then go to next
682  if (roiFirstBinTick > stopTick || roiLastBinTick < startTick) continue;
683 
684  wireRangePtr = &range;
685  break;
686  }
687 
688  // Check that we have found the wire range
689  // Note that if we have not matched an ROI then we can't have a hit either so skip search for that...
690  if (wireRangePtr)
691  {
692  const recob::Hit* rejectedHit = 0;
693  const recob::Hit* bestHit = 0;
694 
695  nMatchedWires++;
696 
697  // The next mission is to recover the hits associated to this Wire
698  // The easiest way to do this is to simply look up all the hits on this channel and then match
699  ChanToHitVecMap::iterator hitIter = channelToHitVec.find(chanToTDCToIDEMap.first);
700 
701  if (hitIter != channelToHitVec.end())
702  {
703  // Loop through the hits for this channel and look for matches
704  // In the event of more than one hit associated to the sim channel range, keep only
705  // the best match (assuming the nearby hits are "extra")
706  // Note that assumption breaks down for long pulse trains but worry about that later
707  for(const auto& hit : hitIter->second)
708  {
709  unsigned short hitStartTick = hit->PeakTime() - fSigmaVec[plane] * hit->RMS();
710  unsigned short hitStopTick = hit->PeakTime() + fSigmaVec[plane] * hit->RMS();
711 
712  // If hit is out of range then skip, it is not related to this particle
713  if (hitStartTick > stopTick || hitStopTick < startTick)
714  {
715  nFakeHitVec[plane]++;
716  rejectedHit = hit;
717  continue;
718  }
719 
720  float hitHeight = hit->PeakAmplitude();
721 
722  // Use the hit with the largest pulse height as the "best"
723  if (hitHeight < hitPeakAmpBest) continue;
724 
725  hitPeakAmpBest = hitHeight;
726  bestHit = hit;
727  hitStartTickBest = hitStartTick;
728  hitStopTickBest = hitStopTick;
729  }
730 
731  // Find a match?
732  if (bestHit)
733  {
734  // std::cout << "found a bestHit (L720)" << std::endl;
735  nElectronsTotalBest = 0.;
736  hitPeakTimeBest = bestHit->PeakTime();
737  hitIntegralBest = bestHit->Integral();
738  hitSummedADCBest = bestHit->SummedADC();
739  hitRMSBest = bestHit->RMS();
740  hitMultiplicityBest = bestHit->Multiplicity();
741  hitLocalIndexBest = bestHit->LocalIndex();
742  hitGoodnessBest = bestHit->GoodnessOfFit();
743  hitNumDegreesBest = bestHit->DegreesOfFreedom();
744  hitSnippetLenBest = bestHit->EndTick() - bestHit->StartTick();
745  hitBaselineBest = 0.; // To do...
746 
747  nMatchedHits++;
748 
749  // Get the number of electrons
750  for(unsigned short tick = hitStartTickBest; tick <= hitStopTickBest; tick++)
751  {
752  unsigned short hitTDC = clockData.TPCTick2TDC(tick - fOffsetVec[plane]);
753 
754  TDCToIDEMap::iterator ideIterator = tdcToIDEMap.find(hitTDC);
755 
756  if (ideIterator != tdcToIDEMap.end()) nElectronsTotalBest += ideIterator->second.numElectrons;
757  }
758  }
759 
760  if (nMatchedHits > 0)
761  {
762  float chgRatioADC = hitSummedADCBest > 0. ? totalElectrons / hitSummedADCBest : 0.;
763  float chgRatioInt = hitIntegralBest > 0. ? totalElectrons / hitIntegralBest : 0.;
764 
765  fHitSumADCVec[plane]->Fill(hitSummedADCBest, 1.);
766  fHitIntegralHistVec[plane]->Fill(hitIntegralBest, 1.);
767  fSimDivHitChgVec[plane]->Fill(chgRatioADC, 1.);
768  fSimDivHitChg1Vec[plane]->Fill(chgRatioInt, 1.);
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.);
772  fHitPulseHeightVec[plane]->Fill(std::min(hitPeakAmpBest,float(149.5)), 1.);
773  fHitPulseWidthVec[plane]->Fill(std::min(hitRMSBest,float(19.8)), 1.);
774  fHitElectronsVec[plane]->Fill(nElectronsTotalBest, 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.);
777  fDeltaMidTDCVec[plane]->Fill(hitPeakTimeBest - maxETick, 1.);
778 
779  nRecobHitVec[plane]++;
780 
781  }
782  else if (rejectedHit)
783  {
784  unsigned short hitStartTick = rejectedHit->PeakTime() - fSigmaVec[plane] * rejectedHit->RMS();
785  unsigned short hitStopTick = rejectedHit->PeakTime() + fSigmaVec[plane] * rejectedHit->RMS();
786 
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;
789  }
790  else
791  {
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;
793  }
794  }
795  }
796  }
797 
798  fWireEfficVec.at(plane)->Fill(totalElectrons, std::min(nMatchedWires,1), 1.);
799  fWireEfficPHVec.at(plane)->Fill(maxElectrons, std::min(nMatchedWires,1), 1.);
800 
801  float matchHit = std::min(nMatchedHits,1);
802  float snippetLen = std::min(float(stopTick - startTick),float(99.5));
803 
804  fNMatchedHitVec[plane]->Fill(nMatchedHits, 1.);
805  fHitEfficVec[plane]->Fill(totalElectrons, matchHit, 1.);
806  fHitEfficPHVec[plane]->Fill(maxElectrons, matchHit, 1.);
807  fHitEfficXZVec[plane]->Fill(cosThetaXZ, matchHit, 1.);
808  fCosXZvRMSVec[plane]->Fill(cosThetaXZ, snippetLen, 1.);
809  fHitEfficRMSVec[plane]->Fill(snippetLen, matchHit, 1.);
810 
811  fHitENEvXZVec[plane]->Fill(cosThetaXZ, totalElectrons, matchHit);
812 
813  // Store tuple variables
814  fTPCVec.push_back(wids[0].TPC);
815  fCryoVec.push_back(wids[0].Cryostat);
816  fPlaneVec.push_back(wids[0].Plane);
817  fWireVec.push_back(wids[0].Wire);
818 
819  fTotalElectronsVec.push_back(totalElectrons);
820  fMaxElectronsVec.push_back(maxElectrons);
821  fStartTickVec.push_back(startTick);
822  fStopTickVec.push_back(stopTick);
823  fMaxETickVec.push_back(maxETick);
824  fPartDirX.push_back(partDirVec[0]);
825  fPartDirY.push_back(partDirVec[1]);
826  fPartDirZ.push_back(partDirVec[2]);
827 
828  fNMatchedWires.push_back(nMatchedWires);
829  fNMatchedHits.push_back(nMatchedHits);
830 
831  fHitPeakTimeVec.push_back(hitPeakTimeBest);
832  fHitPeakAmpVec.push_back(hitPeakAmpBest);
833  fHitPeakRMSVec.push_back(hitRMSBest);
834  fHitBaselinevec.push_back(hitBaselineBest);
835  fHitSummedADCVec.push_back(hitSummedADCBest);
836  fHitIntegralVec.push_back(hitIntegralBest);
837  fHitStartTickVec.push_back(hitStartTickBest);
838  fHitStopTickVec.push_back(hitStopTickBest);
839  fHitMultiplicityVec.push_back(hitMultiplicityBest);
840  fHitLocalIndexVec.push_back(hitLocalIndexBest);
841  fHitGoodnessVec.push_back(hitGoodnessBest);
842  fNumDegreesVec.push_back(hitNumDegreesBest);
843  }
844  }
845 
846  for(size_t idx = 0; idx < fGeometry->Nplanes();idx++)
847  {
848  if (nSimChannelHitVec[idx] > 10)
849  {
850  float hitEfficiency = float(nRecobHitVec[idx]) / float(nSimChannelHitVec[idx]);
851 
852  fNSimChannelHitsVec[idx]->Fill(std::min(nSimChannelHitVec[idx],1999),1.);
853  fNRecobHitVec[idx]->Fill(std::min(nRecobHitVec[idx],1999), 1.);
854  fNFakeHitVec[idx]->Fill(nFakeHitVec[idx]/(float)nSimulatedWiresVec[idx],1.);
855  fHitEfficiencyVec[idx]->Fill(hitEfficiency, 1.);
856  }
857  }
858 
859  return;
860 }
861 
862 // Useful for normalizing histograms
863 void TrackHitEfficiencyAnalysis::endJob(int numEvents)
864 {
865  return;
866 }
867 
868 DEFINE_ART_CLASS_TOOL(TrackHitEfficiencyAnalysis)
869 }
short int LocalIndex() const
How well do we believe we know this hit?
Definition: Hit.h:227
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.
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
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.
void initializeTuple(TTree *) override
Interface for initializing the tuple variables.
tick_as<> tick
Tick number, represented by std::ptrdiff_t.
Definition: electronics.h:75
virtual Status_t Status(raw::ChannelID_t channel) const
Returns a status integer with arbitrary meaning.
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
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
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