All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
HitFinderAnalysis_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 
36 {
37  ////////////////////////////////////////////////////////////////////////
38  //
39  // Class: HitFinderAnalysis
40  // Module Type: producer
41  // File: HitFinderAnalysis.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 HitFinderAnalysis(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  art::InputTag fRawDigitProducerLabel;
111  art::InputTag fWireProducerLabel;
112  art::InputTag fHitProducerLabel;
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*> fHitEfficiencyVec;
152 
153  // TTree variables
154  mutable TTree* fTree;
155 
156  mutable std::vector<int> fTPCVec;
157  mutable std::vector<int> fCryoVec;
158  mutable std::vector<int> fPlaneVec;
159  mutable std::vector<int> fWireVec;
160 
161  mutable std::vector<float> fTotalElectronsVec;
162  mutable std::vector<float> fMaxElectronsVec;
163  mutable std::vector<int> fStartTickVec;
164  mutable std::vector<int> fStopTickVec;
165  mutable std::vector<int> fMaxETickVec;
166  mutable std::vector<float> fPartDirX;
167  mutable std::vector<float> fPartDirY;
168  mutable std::vector<float> fPartDirZ;
169 
170  mutable std::vector<int> fNMatchedWires;
171  mutable std::vector<int> fNMatchedHits;
172 
173  mutable std::vector<float> fHitPeakTimeVec;
174  mutable std::vector<float> fHitPeakAmpVec;
175  mutable std::vector<float> fHitPeakRMSVec;
176  mutable std::vector<float> fHitBaselinevec;
177  mutable std::vector<float> fHitSummedADCVec;
178  mutable std::vector<float> fHitIntegralVec;
179  mutable std::vector<int> fHitStartTickVec;
180  mutable std::vector<int> fHitStopTickVec;
181  mutable std::vector<int> fHitMultiplicityVec;
182  mutable std::vector<int> fHitLocalIndexVec;
183  mutable std::vector<float> fHitGoodnessVec;
184  mutable std::vector<int> fNumDegreesVec;
185 
186  // Useful services, keep copies for now (we can update during begin run periods)
187  const geo::GeometryCore* fGeometry; ///< pointer to Geometry service
188 };
189 
190 //----------------------------------------------------------------------------
191 /// Constructor.
192 ///
193 /// Arguments:
194 ///
195 /// pset - Fcl parameters.
196 ///
197 HitFinderAnalysis::HitFinderAnalysis(fhicl::ParameterSet const & pset) : fTree(nullptr)
198 {
199  fGeometry = lar::providerFrom<geo::Geometry>();
200 
201  configure(pset);
202 
203  // Report.
204  mf::LogInfo("HitFinderAnalysis") << "HitFinderAnalysis configured\n";
205 }
206 
207 //----------------------------------------------------------------------------
208 /// Destructor.
209 HitFinderAnalysis::~HitFinderAnalysis()
210 {}
211 
212 //----------------------------------------------------------------------------
213 /// Reconfigure method.
214 ///
215 /// Arguments:
216 ///
217 /// pset - Fcl parameter set.
218 ///
219 void HitFinderAnalysis::configure(fhicl::ParameterSet const & pset)
220 {
221  fRawDigitProducerLabel = pset.get< art::InputTag >("RawDigitLabel", "rawdigitfilter");
222  fWireProducerLabel = pset.get< art::InputTag >("WireModuleLabel", "decon1droi");
223  fHitProducerLabel = pset.get< art::InputTag >("HitModuleLabel", "gauss");
224  fMCParticleProducerLabel = pset.get< art::InputTag >("MCParticleLabel", "largeant");
225  fSimChannelProducerLabel = pset.get< art::InputTag >("SimChannelLabel", "largeant");
226  fBadChannelProducerLabel = pset.get< art::InputTag >("BadChannelLabel", "simnfspl1:badchannels");
227  fUseBadChannelDB = pset.get< bool >("UseBadChannelDB", true);
228  fLocalDirName = pset.get<std::string >("LocalDirName", std::string("wow"));
229  fOffsetVec = pset.get<std::vector<int> >("OffsetVec", std::vector<int>()={0,0,0});
230  fSigmaVec = pset.get<std::vector<float> >("SigmaVec", std::vector<float>()={1.,1.,1.});
231  fMinAllowedChanStatus = pset.get< int >("MinAllowedChannelStatus");
232 }
233 
234 //----------------------------------------------------------------------------
235 /// Begin job method.
236 void HitFinderAnalysis::initializeHists(art::ServiceHandle<art::TFileService>& tfs, const std::string& dirName)
237 {
238  // Make a directory for these histograms
239  art::TFileDirectory dir = tfs->mkdir(dirName.c_str());
240 
244  fHitSumADCVec.resize(fGeometry->Nplanes());
248  fSimNumTDCVec.resize(fGeometry->Nplanes());
249  fHitNumTDCVec.resize(fGeometry->Nplanes());
250  fSnippetLenVec.resize(fGeometry->Nplanes());
251  fNMatchedHitVec.resize(fGeometry->Nplanes());
252  fDeltaMidTDCVec.resize(fGeometry->Nplanes());
253  fHitVsSimChgVec.resize(fGeometry->Nplanes());
254  fHitVsSimIntVec.resize(fGeometry->Nplanes());
255  fCosXZvRMSVec.resize(fGeometry->Nplanes());
257  fHitENEvXZVec.resize(fGeometry->Nplanes());
259  fNRecobHitVec.resize(fGeometry->Nplanes());
263 
264  fWireEfficVec.resize(fGeometry->Nplanes());
265  fWireEfficPHVec.resize(fGeometry->Nplanes());
266 
267  fHitEfficVec.resize(fGeometry->Nplanes());
268  fHitEfficPHVec.resize(fGeometry->Nplanes());
269  fHitEfficXZVec.resize(fGeometry->Nplanes());
270  fHitEfficRMSVec.resize(fGeometry->Nplanes());
271 
272  for(size_t plane = 0; plane < fGeometry->Nplanes(); plane++)
273  {
274  fTotalElectronsHistVec.at(plane) = dir.make<TH1F>(("TotalElecs" + std::to_string(plane)).c_str(), ";Total # electrons", 250, 0., 100000.);
275  fMaxElectronsHistVec.at(plane) = dir.make<TH1F>(("MaxElecs" + std::to_string(plane)).c_str(), ";Max # electrons", 250, 0., 20000.);
276  fHitElectronsVec.at(plane) = dir.make<TH1F>(("HitElecs" + std::to_string(plane)).c_str(), ";# e- in Hit Range", 250, 0., 100000.);
277  fHitSumADCVec.at(plane) = dir.make<TH1F>(("SumADC" + std::to_string(plane)).c_str(), "Hit Sum ADC", 200, 0., 1000.);
278  fHitIntegralHistVec.at(plane) = dir.make<TH1F>(("Integral" + std::to_string(plane)).c_str(), "Hit Integral", 200, 0., 1000.);
279  fHitPulseHeightVec.at(plane) = dir.make<TH1F>(("PulseHeight" + std::to_string(plane)).c_str(), "Hit PH (ADC)", 200, 0., 150.);
280  fHitPulseWidthVec.at(plane) = dir.make<TH1F>(("PulseWidth" + std::to_string(plane)).c_str(), "Hit PulseWidth;Hit RMS", 100, 0., 20.);
281  fSimNumTDCVec.at(plane) = dir.make<TH1F>(("SimNumTDC" + std::to_string(plane)).c_str(), ";TDC ticks", 100, 0., 100.);
282  fHitNumTDCVec.at(plane) = dir.make<TH1F>(("HitNumTDC" + std::to_string(plane)).c_str(), ";ticks", 100, 0., 100.);
283  fSnippetLenVec.at(plane) = dir.make<TH1F>(("SnippetLen" + std::to_string(plane)).c_str(), ";ticks", 100, 0., 100.);
284  fNMatchedHitVec.at(plane) = dir.make<TH1F>(("NMatched" + std::to_string(plane)).c_str(), ";# hits", 20, 0., 20.);
285  fDeltaMidTDCVec.at(plane) = dir.make<TH1F>(("DeltaMid" + std::to_string(plane)).c_str(), ";# hits", 50, -25., 25.);
286  fNSimChannelHitsVec.at(plane) = dir.make<TH1F>(("NSimChan" + std::to_string(plane)).c_str(), ";# hits", 100, 0., 1200.);
287  fNRecobHitVec.at(plane) = dir.make<TH1F>(("NRecobHit" + std::to_string(plane)).c_str(), ";# hits", 100, 0., 1200.);
288  fHitEfficiencyVec.at(plane) = dir.make<TH1F>(("PlnEffic" + std::to_string(plane)).c_str(), ";# hits", 101, 0., 1.01);
289  fSimDivHitChgVec.at(plane) = dir.make<TH1F>(("SimDivHit" + std::to_string(plane)).c_str(), ";# e / SummedADC", 200, 0., 200.);
290  fSimDivHitChg1Vec.at(plane) = dir.make<TH1F>(("SimDivHit1" + std::to_string(plane)).c_str(), ";# e / Integral", 200, 0., 200.);
291 
292  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.);
293  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.);
294  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.);
295 
296  fWireEfficVec.at(plane) = dir.make<TProfile> (("WireEffic" + std::to_string(plane)).c_str(), "Wire Efficiency;# electrons", 50, 0., 100000., 0., 1.);
297  fWireEfficPHVec.at(plane) = dir.make<TProfile> (("WireEfficPH" + std::to_string(plane)).c_str(), "Wire Efficiency;# electrons", 50, 0., 20000., 0., 1.);
298 
299  fHitEfficVec.at(plane) = dir.make<TProfile> (("HitEffic" + std::to_string(plane)).c_str(), "Hit Efficiency;Total # e-", 50, 0., 100000., 0., 1.);
300  fHitEfficPHVec.at(plane) = dir.make<TProfile> (("HitEfficPH" + std::to_string(plane)).c_str(), "Hit Efficiency;Max # e-", 50, 0., 20000., 0., 1.);
301  fHitEfficXZVec.at(plane) = dir.make<TProfile> (("HitEfficXZ" + std::to_string(plane)).c_str(), "Hit Efficiency;cos(thetaXZ)", 50, -1., 1., 0., 1.);
302  fHitEfficRMSVec.at(plane) = dir.make<TProfile> (("HitEfficRMS" + std::to_string(plane)).c_str(), "Hit Efficiency;MC Length (ticks)", 50, 0., 100., 0., 1.);
303  fCosXZvRMSVec.at(plane) = dir.make<TProfile> (("CosXZVRMS" + std::to_string(plane)).c_str(), "CosXZ v. RMS;Cos(XZ);Ticks", 50, -1., 1., 0., 100.);
304  fHitENEvXZVec.at(plane) = dir.make<TProfile2D>(("HitENEvXZ" + std::to_string(plane)).c_str(), "XZ v # e;cos(thetaXZ);#electrons", 50, -1., 1.,
305  50, 0., 100000., 0., 1.);
306  }
307 
308  return;
309 }
310 
311 void HitFinderAnalysis::initializeTuple(TTree* tree)
312 {
313  fTree = tree;
314 
315  fTree->Branch("CryostataVec", "std::vector<int>", &fCryoVec);
316  fTree->Branch("TPCVec", "std::vector<int>", &fTPCVec);
317  fTree->Branch("PlaneVec", "std::vector<int>", &fPlaneVec);
318  fTree->Branch("WireVec", "std::vector<int>", &fWireVec);
319 
320  fTree->Branch("TotalElectronsVec", "std::vector<float>", &fTotalElectronsVec);
321  fTree->Branch("MaxElectronsVec", "std::vector<float>", &fMaxElectronsVec);
322  fTree->Branch("StartTick", "std::vector<int>", &fStartTickVec);
323  fTree->Branch("StopTick", "std::vector<int>", &fStopTickVec);
324  fTree->Branch("MaxETick", "std::vector<int>", &fMaxETickVec);
325  fTree->Branch("PartDirX", "std::vector<float>", &fPartDirX);
326  fTree->Branch("PartDirY", "std::vector<float>", &fPartDirY);
327  fTree->Branch("PartDirZ", "std::vector<float>", &fPartDirZ);
328 
329  fTree->Branch("NMatchedWires", "std::vector<int>", &fNMatchedWires);
330  fTree->Branch("NMatchedHits", "std::vector<int>", &fNMatchedHits);
331 
332  fTree->Branch("HitPeakTimeVec", "std::vector<float>", &fHitPeakTimeVec);
333  fTree->Branch("HitPeakAmpVec", "std::vector<float>", &fHitPeakAmpVec);
334  fTree->Branch("HitPeakRMSVec", "std::vector<float>", &fHitPeakRMSVec);
335  fTree->Branch("HitBaselineVec", "std::vector<float>", &fHitBaselinevec);
336  fTree->Branch("HitSummedADCVec", "std::vector<float>", &fHitSummedADCVec);
337  fTree->Branch("HitIntegralVec", "std::vector<float>", &fHitIntegralVec);
338  fTree->Branch("HitStartTickVec", "std::vector<int>", &fHitStartTickVec);
339  fTree->Branch("HitStopTickVec", "std::vector<int>", &fHitStopTickVec);
340  fTree->Branch("HitMultiplicity", "std::vector<int>", &fHitMultiplicityVec);
341  fTree->Branch("HitLocalIndex", "std::vector<int>", &fHitLocalIndexVec);
342  fTree->Branch("HitGoodness", "std::vector<float>", &fHitGoodnessVec);
343  fTree->Branch("HitNumDegrees", "std::vector<int>", &fNumDegreesVec);
344 
345  clear();
346 
347  return;
348 }
349 
350 void HitFinderAnalysis::clear() const
351 {
352  fTPCVec.clear();
353  fCryoVec.clear();
354  fPlaneVec.clear();
355  fWireVec.clear();
356 
357  fTotalElectronsVec.clear();
358  fMaxElectronsVec.clear();
359  fStartTickVec.clear();
360  fStopTickVec.clear();
361  fMaxETickVec.clear();
362  fPartDirX.clear();
363  fPartDirY.clear();
364  fPartDirZ.clear();
365 
366  fNMatchedWires.clear();
367  fNMatchedHits.clear();
368 
369  fHitPeakTimeVec.clear();
370  fHitPeakAmpVec.clear();
371  fHitPeakRMSVec.clear();
372  fHitBaselinevec.clear();
373  fHitSummedADCVec.clear();
374  fHitIntegralVec.clear();
375  fHitStartTickVec.clear();
376  fHitStopTickVec.clear();
377  fHitMultiplicityVec.clear();
378  fHitLocalIndexVec.clear();
379  fHitGoodnessVec.clear();
380  fNumDegreesVec.clear();
381 
382  return;
383 }
384 
385 void HitFinderAnalysis::fillHistograms(const art::Event& event) const
386 {
387  // Always clear the tuple
388  clear();
389 
390  art::Handle< std::vector<sim::SimChannel>> simChannelHandle;
391  event.getByLabel(fSimChannelProducerLabel, simChannelHandle);
392 
393  // If there is no sim channel informaton then exit
394  if (!simChannelHandle.isValid() || simChannelHandle->empty()) return;
395 
396  art::Handle< std::vector<raw::RawDigit> > rawDigitHandle;
397  event.getByLabel(fRawDigitProducerLabel, rawDigitHandle);
398 
399  art::Handle< std::vector<recob::Wire> > wireHandle;
400  event.getByLabel(fWireProducerLabel, wireHandle);
401 
402  art::Handle< std::vector<recob::Hit> > hitHandle;
403  event.getByLabel(fHitProducerLabel, hitHandle);
404 
405  art::Handle< std::vector<simb::MCParticle>> mcParticleHandle;
406  event.getByLabel(fMCParticleProducerLabel, mcParticleHandle);
407 
408  if (!rawDigitHandle.isValid() || !wireHandle.isValid() || !hitHandle.isValid() || !mcParticleHandle.isValid()) return;
409 
410  // Look up the list of bad channels
411  art::Handle< std::vector<int>> badChannelHandle;
412  event.getByLabel(fBadChannelProducerLabel, badChannelHandle);
413 
414  // Find the associations between wire data and hits
415  // What we want to be able to do is look up hits that have been associated to Wire data
416  // So we ask for the list of hits, given a handle to the Wire data and remember that the
417  // associations are made in the hit finder
418  //art::FindManyP<recob::Hit> wireHitAssns(wireHandle,event,fHitProducerLabel);
419 
420  // what needs to be done?
421  // First we should build out a straightforward channel to Wire map so we can look up a given
422  // channel's Wire data as we loop over SimChannels.
423  using ChanToWireMap = std::unordered_map<raw::ChannelID_t,const recob::Wire*>;
424 
425  ChanToWireMap channelToWireMap;
426 
427  for(const auto& wire : *wireHandle) channelToWireMap[wire.Channel()] = &wire;
428 
429  // We will use the presence of a RawDigit as an indicator of a good channel... So
430  // we want a mapping between channel and RawDigit
431  using ChanToRawDigitMap = std::unordered_map<raw::ChannelID_t,const raw::RawDigit*>;
432 
433  ChanToRawDigitMap chanToRawDigitMap;
434 
435  for(const auto& rawDigit : *rawDigitHandle) chanToRawDigitMap[rawDigit.Channel()] = &rawDigit;
436 
437  // Then we should map out all hits by channel so we can easily look up from sim channels
438  // Then go through the sim channels and match hits
439  using ChanToHitVecMap = std::unordered_map<raw::ChannelID_t,std::vector<const recob::Hit*>>;
440  ChanToHitVecMap channelToHitVec;
441 
442  for(const auto& hit : *hitHandle) channelToHitVec[hit.Channel()].push_back(&hit);
443 
444  // It is useful to create a mapping between trackID and MCParticle
445  using TrackIDToMCParticleMap = std::unordered_map<int, const simb::MCParticle*>;
446 
447  TrackIDToMCParticleMap trackIDToMCParticleMap;
448 
449  for(const auto& mcParticle : *mcParticleHandle) trackIDToMCParticleMap[mcParticle.TrackId()] = &mcParticle;
450 
451  // Now go through the sim channels
452  // There are several things going on here... for each channel we have particles (track id's) depositing energy in a range to ticks
453  // So... for each channel we want to build a structure that relates particles to tdc ranges and deposited energy (or electrons)
454  // Here is a complicated structure:
455 
456  using TDCToIDEMap = std::map<unsigned short, sim::IDE>; // We need this one in order
457  using ChanToTDCToIDEMap = std::map<raw::ChannelID_t, TDCToIDEMap>;
458  using PartToChanToTDCToIDEMap = std::unordered_map<int, ChanToTDCToIDEMap>;
459 
460  PartToChanToTDCToIDEMap partToChanToTDCToIDEMap;
461 
462  // Build out the above data structure
463  for(const auto& simChannel : *simChannelHandle)
464  {
465  for(const auto& tdcide : simChannel.TDCIDEMap())
466  {
467  for(const auto& ide : tdcide.second) partToChanToTDCToIDEMap[ide.trackID][simChannel.Channel()][tdcide.first] = ide;
468  }
469  }
470 
471  auto const clockData = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(event);
472  const lariov::ChannelStatusProvider& chanFilt = art::ServiceHandle<lariov::ChannelStatusService>()->GetProvider();
473 
474  std::vector<int> nSimChannelHitVec = {0,0,0};
475  std::vector<int> nRecobHitVec = {0,0,0};
476 
477  for(const auto& partToChanInfo : partToChanToTDCToIDEMap)
478  {
479  TrackIDToMCParticleMap::const_iterator trackIDToMCPartItr = trackIDToMCParticleMap.find(partToChanInfo.first);
480 
481  if (trackIDToMCPartItr == trackIDToMCParticleMap.end()) continue;
482 
483  int trackPDGCode = trackIDToMCPartItr->second->PdgCode();
484  std::string processName = trackIDToMCPartItr->second->Process();
485 
486  // Looking for primary muons (e.g. CR Tracks)
487  if (fabs(trackPDGCode) != 13 || processName != "primary") continue;
488 
489  // Recover particle position and angle information
490  Eigen::Vector3f partStartPos(trackIDToMCPartItr->second->Vx(),trackIDToMCPartItr->second->Vy(),trackIDToMCPartItr->second->Vz());
491  Eigen::Vector3f partStartDir(trackIDToMCPartItr->second->Px(),trackIDToMCPartItr->second->Py(),trackIDToMCPartItr->second->Pz());
492 
493  partStartDir.normalize();
494 
495  Eigen::Vector2f partStartDirVecXZ(partStartDir[0],partStartDir[2]);
496 
497  partStartDirVecXZ.normalize();
498 
499  // Assuming the SimChannels contain position information (currently not true for WC produced SimChannels)
500  // then we want to keep a running position
501  std::vector<Eigen::Vector3f> lastPositionVec = {partStartPos,partStartPos,partStartPos};
502 
503  for(const auto& chanToTDCToIDEMap : partToChanInfo.second)
504  {
505  // skip bad channels
506  if (fUseBadChannelDB)
507  {
508  // This is the "correct" way to check and remove bad channels...
509  if( chanFilt.Status(chanToTDCToIDEMap.first) < fMinAllowedChanStatus)
510  {
511  //std::vector<geo::WireID> wids = fGeometry->ChannelToWire(chanToTDCToIDEMap.first);
512  //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;
513  continue;
514  }
515  }
516 
517  // Was a list made available?
518  // If so then we try that
519  if (badChannelHandle.isValid())
520  {
521  // Here we query the input list from the wirecell processing
522  std::vector<int>::const_iterator badItr = std::find(badChannelHandle->begin(),badChannelHandle->end(),chanToTDCToIDEMap.first);
523 
524  if (badItr != badChannelHandle->end()) continue;
525  // {
526  // ChanToRawDigitMap::const_iterator rawDigitItr = chanToRawDigitMap.find(chanToTDCToIDEMap.first);
527  //
528  // if (rawDigitItr != chanToRawDigitMap.end())
529  // {
530  // float nSig(3.);
531  // float mean(0.);
532  // float rmsFull(0.);
533  // float rmsTrunc(0.);
534  // int nTrunc(0);
535  //
536  // getTruncatedMeanRMS(rawDigitItr->second->ADCs(), nSig, mean, rmsFull, rmsTrunc, nTrunc);
537  //
538  // std::cout << "--> Rejecting channel: " << chanToTDCToIDEMap.first << " from bad channel list, rms: " << rmsFull << std::endl;
539  // }
540  //
541  // continue;
542  // }
543  }
544 
545  TDCToIDEMap tdcToIDEMap = chanToTDCToIDEMap.second;
546  float totalElectrons(0.);
547  float maxElectrons(0.);
548  unsigned short maxElectronsTDC(0);
549  int nMatchedWires(0);
550  int nMatchedHits(0);
551 
552  // The below try-catch block may no longer be necessary
553  // Decode the channel and make sure we have a valid one
554  std::vector<geo::WireID> wids = fGeometry->ChannelToWire(chanToTDCToIDEMap.first);
555 
556  // Recover plane and wire in the plane
557  unsigned int plane = wids[0].Plane;
558 // unsigned int wire = wids[0].Wire;
559 
560  Eigen::Vector3f avePosition(0.,0.,0.);
561 
562  for(const auto& ideVal : tdcToIDEMap)
563  {
564  totalElectrons += ideVal.second.numElectrons;
565 
566  if (maxElectrons < ideVal.second.numElectrons)
567  {
568  maxElectrons = ideVal.second.numElectrons;
569  maxElectronsTDC = ideVal.first;
570  }
571 
572  avePosition += Eigen::Vector3f(ideVal.second.x,ideVal.second.y,ideVal.second.z);
573  }
574 
575  // Get local track direction by using the average position of deposited charge as the current position
576  // and then subtracting the last position
577  avePosition /= float(tdcToIDEMap.size());
578 
579  Eigen::Vector3f partDirVec = avePosition - lastPositionVec[plane];
580 
581  partDirVec.normalize();
582 
583  lastPositionVec[plane] = avePosition;
584 
585  // Now what we want is the projected angle in the XZ plane
586  Eigen::Vector2f projPairDirVec(partDirVec[0],partDirVec[2]);
587 
588  projPairDirVec.normalize();
589 
590  float cosThetaXZ = projPairDirVec[0];
591 
592  totalElectrons = std::min(totalElectrons, float(99900.));
593 
594  fTotalElectronsHistVec[plane]->Fill(totalElectrons, 1.);
595  fMaxElectronsHistVec[plane]->Fill(maxElectrons, 1.);
596 
597  nSimChannelHitVec[plane]++;
598 
599  unsigned short startTDC = tdcToIDEMap.begin()->first;
600  unsigned short stopTDC = tdcToIDEMap.rbegin()->first;
601 
602  // Convert to ticks to get in same units as hits
603  unsigned short startTick = clockData.TPCTDC2Tick(startTDC) + fOffsetVec[plane];
604  unsigned short stopTick = clockData.TPCTDC2Tick(stopTDC) + fOffsetVec[plane];
605  unsigned short maxETick = clockData.TPCTDC2Tick(maxElectronsTDC) + fOffsetVec[plane];
606 
607  fSimNumTDCVec[plane]->Fill(stopTick - startTick, 1.);
608 
609  // Set up to extract the "best" parameters in the event of more than one hit for this pulse train
610  float nElectronsTotalBest(0.);
611  float hitSummedADCBest(0.);
612  float hitIntegralBest(0.);
613  float hitPeakTimeBest(0.);
614  float hitPeakAmpBest(-100.);
615  float hitRMSBest(0.);
616  int hitMultiplicityBest(0);
617  int hitLocalIndexBest(0);
618  float hitGoodnessBest(0.);
619  int hitNumDegreesBest(0);
620  float hitBaselineBest(0.);
621  float hitSnippetLenBest(0.);
622  unsigned short hitStopTickBest(0);
623  unsigned short hitStartTickBest(0);
624 // unsigned short midHitTickBest(0);
625 
626  // Start by recovering the Wire associated to this channel
627  ChanToWireMap::const_iterator wireItr = channelToWireMap.find(chanToTDCToIDEMap.first);
628 
629  if (wireItr != channelToWireMap.end())
630  {
631  const recob::Wire::RegionsOfInterest_t& signalROI = wireItr->second->SignalROI();
632  const lar::sparse_vector<float>::datarange_t* wireRangePtr(NULL);
633 
634  // Here we need to match the range of the ROI's on the given Wire with the tick range from the SimChannel
635  for(const auto& range : signalROI.get_ranges())
636  {
637  // #################################################
638  // ### Getting a vector of signals for this wire ###
639  // #################################################
640  //std::vector<float> signal(wire->Signal());
641 
642  const std::vector<float>& signal = range.data();
643  raw::TDCtick_t roiFirstBinTick = range.begin_index();
644  raw::TDCtick_t roiLastBinTick = roiFirstBinTick + signal.size();
645 
646  // If no overlap then go to next
647  if (roiFirstBinTick > stopTick || roiLastBinTick < startTick) continue;
648 
649  wireRangePtr = &range;
650  break;
651  }
652 
653  // Check that we have found the wire range
654  // Note that if we have not matched an ROI then we can't have a hit either so skip search for that...
655  if (wireRangePtr)
656  {
657  const recob::Hit* rejectedHit = 0;
658  const recob::Hit* bestHit = 0;
659 
660  nMatchedWires++;
661 
662  // The next mission is to recover the hits associated to this Wire
663  // The easiest way to do this is to simply look up all the hits on this channel and then match
664  ChanToHitVecMap::iterator hitIter = channelToHitVec.find(chanToTDCToIDEMap.first);
665 
666  if (hitIter != channelToHitVec.end())
667  {
668  // Loop through the hits for this channel and look for matches
669  // In the event of more than one hit associated to the sim channel range, keep only
670  // the best match (assuming the nearby hits are "extra")
671  // Note that assumption breaks down for long pulse trains but worry about that later
672  for(const auto& hit : hitIter->second)
673  {
674  unsigned short hitStartTick = hit->PeakTime() - fSigmaVec[plane] * hit->RMS();
675  unsigned short hitStopTick = hit->PeakTime() + fSigmaVec[plane] * hit->RMS();
676 // unsigned short midHitTick = (hitStopTick + hitStartTick) / 2;
677 
678  // If hit is out of range then skip, it is not related to this particle
679  if (hitStartTick > stopTick || hitStopTick < startTick)
680  {
681  if (plane == 2) rejectedHit = hit;
682  continue;
683  }
684 
685  float hitHeight = hit->PeakAmplitude();
686 
687  // Use the hit with the largest pulse height as the "best"
688  if (hitHeight < hitPeakAmpBest) continue;
689 
690  hitPeakAmpBest = hitHeight;
691  bestHit = hit;
692  hitStartTickBest = hitStartTick;
693  hitStopTickBest = hitStopTick;
694 // midHitTickBest = midHitTick;
695  }
696 
697  // Find a match?
698  if (bestHit)
699  {
700  nElectronsTotalBest = 0.;
701  hitPeakTimeBest = bestHit->PeakTime();
702  hitIntegralBest = bestHit->Integral();
703  hitSummedADCBest = bestHit->SummedADC();
704  hitRMSBest = bestHit->RMS();
705  hitMultiplicityBest = bestHit->Multiplicity();
706  hitLocalIndexBest = bestHit->LocalIndex();
707  hitGoodnessBest = bestHit->GoodnessOfFit();
708  hitNumDegreesBest = bestHit->DegreesOfFreedom();
709  hitSnippetLenBest = bestHit->EndTick() - bestHit->StartTick();
710  hitBaselineBest = 0.; // To do...
711 
712  nMatchedHits++;
713 
714  // Get the number of electrons
715  for(unsigned short tick = hitStartTickBest; tick <= hitStopTickBest; tick++)
716  {
717  unsigned short hitTDC = clockData.TPCTick2TDC(tick - fOffsetVec[plane]);
718 
719  TDCToIDEMap::iterator ideIterator = tdcToIDEMap.find(hitTDC);
720 
721  if (ideIterator != tdcToIDEMap.end()) nElectronsTotalBest += ideIterator->second.numElectrons;
722  }
723  }
724 
725  if (nMatchedHits > 0)
726  {
727  float chgRatioADC = hitSummedADCBest > 0. ? totalElectrons / hitSummedADCBest : 0.;
728  float chgRatioInt = hitIntegralBest > 0. ? totalElectrons / hitIntegralBest : 0.;
729 
730  fHitSumADCVec[plane]->Fill(hitSummedADCBest, 1.);
731  fHitIntegralHistVec[plane]->Fill(hitIntegralBest, 1.);
732  fSimDivHitChgVec[plane]->Fill(chgRatioADC, 1.);
733  fSimDivHitChg1Vec[plane]->Fill(chgRatioInt, 1.);
734  fHitVsSimChgVec[plane]->Fill(std::min(hitSummedADCBest,float(999.)), totalElectrons, 1.);
735  fHitVsSimIntVec[plane]->Fill(std::min(hitIntegralBest,float(999.)), totalElectrons, 1.);
736  fToteVHitEIntVec[plane]->Fill(std::min(totalElectrons,float(99999.)),std::min(nElectronsTotalBest,float(99999.)),1.);
737  fHitPulseHeightVec[plane]->Fill(std::min(hitPeakAmpBest,float(149.5)), 1.);
738  fHitPulseWidthVec[plane]->Fill(std::min(hitRMSBest,float(19.8)), 1.);
739  fHitElectronsVec[plane]->Fill(nElectronsTotalBest, 1.);
740  fHitNumTDCVec[plane]->Fill(std::min(float(hitStopTickBest - hitStartTickBest),float(99.5)), 1.);
741  fSnippetLenVec[plane]->Fill(std::min(hitSnippetLenBest, float(99.5)), 1.);
742  fDeltaMidTDCVec[plane]->Fill(hitPeakTimeBest - maxETick, 1.);
743 
744  nRecobHitVec[plane]++;
745  }
746  else if (rejectedHit)
747  {
748  unsigned short hitStartTick = rejectedHit->PeakTime() - fSigmaVec[plane] * rejectedHit->RMS();
749  unsigned short hitStopTick = rejectedHit->PeakTime() + fSigmaVec[plane] * rejectedHit->RMS();
750 
751  mf::LogDebug("HitFinderAnalysis") << "**> TPC: " << rejectedHit->WireID().TPC << ", Plane " << rejectedHit->WireID().Plane << ", wire: " << rejectedHit->WireID().Wire << ", hit start/stop tick: " << hitStartTick << "/" << hitStopTick << ", start/stop ticks: " << startTick << "/" << stopTick << std::endl;
752  mf::LogDebug("HitFinderAnalysis") << " 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;
753  }
754  else
755  {
756  mf::LogDebug("HitFinderAnalysis") << "==> No match, TPC/Plane/Wire: " << "/" << wids[0].TPC << "/" << wids[0].Plane << "/" << wids[0].Wire << ", # electrons: " << totalElectrons << ", startTick: " << startTick << ", stopTick: " << stopTick << std::endl;
757  }
758  }
759  }
760  }
761 
762  fWireEfficVec.at(plane)->Fill(totalElectrons, std::min(nMatchedWires,1), 1.);
763  fWireEfficPHVec.at(plane)->Fill(maxElectrons, std::min(nMatchedWires,1), 1.);
764 
765  float matchHit = std::min(nMatchedHits,1);
766  float snippetLen = std::min(float(stopTick - startTick),float(99.5));
767 
768  fNMatchedHitVec[plane]->Fill(nMatchedHits, 1.);
769  fHitEfficVec[plane]->Fill(totalElectrons, matchHit, 1.);
770  fHitEfficPHVec[plane]->Fill(maxElectrons, matchHit, 1.);
771  fHitEfficXZVec[plane]->Fill(cosThetaXZ, matchHit, 1.);
772  fCosXZvRMSVec[plane]->Fill(cosThetaXZ, snippetLen, 1.);
773  fHitEfficRMSVec[plane]->Fill(snippetLen, matchHit, 1.);
774 
775  fHitENEvXZVec[plane]->Fill(cosThetaXZ, totalElectrons, matchHit);
776 
777  // Store tuple variables
778  fTPCVec.push_back(wids[0].TPC);
779  fCryoVec.push_back(wids[0].Cryostat);
780  fPlaneVec.push_back(wids[0].Plane);
781  fWireVec.push_back(wids[0].Wire);
782 
783  fTotalElectronsVec.push_back(totalElectrons);
784  fMaxElectronsVec.push_back(maxElectrons);
785  fStartTickVec.push_back(startTick);
786  fStopTickVec.push_back(stopTick);
787  fMaxETickVec.push_back(maxETick);
788  fPartDirX.push_back(partDirVec[0]);
789  fPartDirY.push_back(partDirVec[1]);
790  fPartDirZ.push_back(partDirVec[2]);
791 
792  fNMatchedWires.push_back(nMatchedWires);
793  fNMatchedHits.push_back(nMatchedHits);
794 
795  fHitPeakTimeVec.push_back(hitPeakTimeBest);
796  fHitPeakAmpVec.push_back(hitPeakAmpBest);
797  fHitPeakRMSVec.push_back(hitRMSBest);
798  fHitBaselinevec.push_back(hitBaselineBest);
799  fHitSummedADCVec.push_back(hitSummedADCBest);
800  fHitIntegralVec.push_back(hitIntegralBest);
801  fHitStartTickVec.push_back(hitStartTickBest);
802  fHitStopTickVec.push_back(hitStopTickBest);
803  fHitMultiplicityVec.push_back(hitMultiplicityBest);
804  fHitLocalIndexVec.push_back(hitLocalIndexBest);
805  fHitGoodnessVec.push_back(hitGoodnessBest);
806  fNumDegreesVec.push_back(hitNumDegreesBest);
807  }
808  }
809 
810  for(size_t idx = 0; idx < fGeometry->Nplanes();idx++)
811  {
812  if (nSimChannelHitVec[idx] > 10)
813  {
814  float hitEfficiency = float(nRecobHitVec[idx]) / float(nSimChannelHitVec[idx]);
815 
816  fNSimChannelHitsVec[idx]->Fill(std::min(nSimChannelHitVec[idx],1999),1.);
817  fNRecobHitVec[idx]->Fill(std::min(nRecobHitVec[idx],1999), 1.);
818  fHitEfficiencyVec[idx]->Fill(hitEfficiency, 1.);
819  }
820  }
821 
822  return;
823 }
824 
825 // Useful for normalizing histograms
826 void HitFinderAnalysis::endJob(int numEvents)
827 {
828  return;
829 }
830 
831 DEFINE_ART_CLASS_TOOL(HitFinderAnalysis)
832 }
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
void configure(fhicl::ParameterSet const &pset) override
std::vector< int > fOffsetVec
Allow offsets for each 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.
void endJob(int numEvents) override
Interface for method to executve at the end of run processing.
std::vector< float > fSigmaVec
Window size for matching to SimChannels.
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
void fillHistograms(const art::Event &) const override
Interface for filling histograms.
float PeakAmplitude() const
The estimated amplitude of the hit at its peak, in ADC units.
Definition: Hit.h:221
int TDCtick_t
Type representing a TDC tick.
Definition: RawTypes.h:25
unsigned int Nplanes(unsigned int tpc=0, unsigned int cstat=0) const
Returns the total number of wire planes in the specified TPC.
HitFinderAnalysis(fhicl::ParameterSet const &pset)
Constructor.
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.
std::string fLocalDirName
Fraction for truncated mean.
const geo::GeometryCore * fGeometry
pointer to Geometry service
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
std::vector< art::Ptr< recob::Hit >> HitPtrVec
tuple dir
Definition: dropbox.py:28
raw::TDCtick_t EndTick() const
Final tdc tick for hit.
Definition: Hit.h:217
std::map< size_t, HitPtrVec > ViewHitMap
void initializeHists(art::ServiceHandle< art::TFileService > &, const std::string &) override
Interface for initializing the histograms to be filled.
float PeakTime() const
Time of the signal peak, in tick units.
Definition: Hit.h:218
std::map< int, ViewHitMap > TrackViewHitMap
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
int fMinAllowedChanStatus
Don&#39;t consider channels with lower status.
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
void initializeTuple(TTree *) override
Interface for initializing the tuple variables.