All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
BasicHitAnalysis_tool.cc
Go to the documentation of this file.
1 #include "icaruscode/Analysis/tools/IHitHistogramTool.h"
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 
11 #include "larcore/Geometry/Geometry.h"
12 #include "larcore/CoreUtils/ServiceUtil.h" // lar::providerFrom()
13 
14 #include "lardataobj/RecoBase/Hit.h"
15 
16 #include "TH1.h"
17 #include "TH2.h"
18 #include "TProfile.h"
19 #include "TProfile2D.h"
20 
21 #include <cmath>
22 #include <algorithm>
23 
24 namespace BasicHitAnalysis
25 {
26  ////////////////////////////////////////////////////////////////////////
27  //
28  // Class: BasicHitAnalysis
29  // Module Type: producer
30  // File: BasicHitAnalysis.h
31  //
32  // The intent of this module is to provide methods for
33  // "analyzing" hits on waveforms
34  //
35  // Configuration parameters:
36  //
37  // TruncMeanFraction - the fraction of waveform bins to discard when
38  //
39  // Created by Tracy Usher (usher@slac.stanford.edu) on February 19, 2016
40  //
41  ////////////////////////////////////////////////////////////////////////
42 
43 // The following typedefs will, obviously, be useful
44 using HitPtrVec = std::vector<art::Ptr<recob::Hit>>;
45 using ViewHitMap = std::map<size_t,HitPtrVec>;
46 using TrackViewHitMap = std::map<int,ViewHitMap>;
47 
48 class BasicHitAnalysis : virtual public IHitHistogramTool
49 {
50 public:
51  /**
52  * @brief Constructor
53  *
54  * @param pset
55  */
56  explicit BasicHitAnalysis(fhicl::ParameterSet const & pset);
57 
58  /**
59  * @brief Destructor
60  */
61  ~BasicHitAnalysis();
62 
63  // provide for initialization
64  void configure(fhicl::ParameterSet const & pset) override;
65 
66  /**
67  * @brief Interface for initializing the histograms to be filled
68  *
69  * @param TFileService handle to the TFile service
70  * @param string subdirectory to store the hists in
71  */
72  void initializeHists(art::ServiceHandle<art::TFileService>&, const std::string&) override;
73 
74  /**
75  * @brief Interface for method to executve at the end of run processing
76  *
77  * @param int number of events to use for normalization
78  */
79  void endJob(int numEvents) override;
80 
81  /**
82  * @brief Interface for filling histograms
83  */
84  void fillHistograms(const HitPtrVec&) const override;
85 
86 private:
87 
88  // Fcl parameters.
89  std::string fLocalDirName; ///< Fraction for truncated mean
90 
91  // Pointers to the histograms we'll create.
92  std::vector<TH1D*> fHitsByWire;
93  TH1D* fDriftTimes[3];
94  TH1D* fHitsByTime[3];
95  TH1D* fPulseHeight[3];
96  TH1D* fPulseHeightSingle[3];
97  TH1D* fPulseHeightMulti[3];
98  TH1D* fPulseHeightMultiMax[3];
99  TH1D* fPulseHeightLong[3];
100  TH1D* fChi2DOF[3];
101  TH1D* fNumDegFree[3];
102  TH1D* fChi2DOFSingle[3];
103  TH1D* fHitMult[3];
104  TH1D* fHitCharge[3];
105  TH1D* fHitChargeSingle[3];
106  TH1D* fHitChargeMulti[3];
107  TH1D* fHitChargeLong[3];
108  TH1D* fFitWidth[3];
109  TH1D* fFitWidthSingle[3];
110  TH1D* fFitWidthMulti[3];
111  TH1D* fFitWidthLong[3];
112  TH1D* fHitSumADC[3];
113  TH2D* fNDFVsChi2[3];
114  TH2D* fPulseHVsWidth[3];
115  TH2D* fPulseHVsCharge[3];
116  TH1D* fBadWPulseHeight;
117  TH2D* fBadWPulseHVsWidth;
118  TH1D* fBadWHitsByWire;
119  TProfile* fPulseHVsHitNo[3];
120  TProfile* fChargeVsHitNo[3];
121  TProfile* fChargeVsHitNoS[3];
122 
123  TH2D* fSPHvsIdx[3];
124  TH2D* fSWidVsIdx[3];
125  TH2D* f1PPHvsWid[3];
126  TH2D* fSPPHvsWid[3];
127  TH2D* fSOPHvsWid[3];
128  TH2D* fPHRatVsIdx[3];
129 
130  // Useful services, keep copies for now (we can update during begin run periods)
131  const geo::GeometryCore* fGeometry; ///< pointer to Geometry service
132 };
133 
134 //----------------------------------------------------------------------------
135 /// Constructor.
136 ///
137 /// Arguments:
138 ///
139 /// pset - Fcl parameters.
140 ///
141 BasicHitAnalysis::BasicHitAnalysis(fhicl::ParameterSet const & pset)
142 {
143  fGeometry = lar::providerFrom<geo::Geometry>();
144 
145  configure(pset);
146 
147  // Report.
148  mf::LogInfo("BasicHitAnalysis") << "BasicHitAnalysis configured\n";
149 }
150 
151 //----------------------------------------------------------------------------
152 /// Destructor.
153 BasicHitAnalysis::~BasicHitAnalysis()
154 {}
155 
156 //----------------------------------------------------------------------------
157 /// Reconfigure method.
158 ///
159 /// Arguments:
160 ///
161 /// pset - Fcl parameter set.
162 ///
163 void BasicHitAnalysis::configure(fhicl::ParameterSet const & pset)
164 {
165  fLocalDirName = pset.get<std::string>("LocalDirName", std::string("wow"));
166 }
167 
168 //----------------------------------------------------------------------------
169 /// Begin job method.
170 void BasicHitAnalysis::initializeHists(art::ServiceHandle<art::TFileService>& tfs, const std::string& dirName)
171 {
172  // Make a directory for these histograms
173  art::TFileDirectory dir = tfs->mkdir(dirName.c_str());
174 
175  fHitsByWire.resize(fGeometry->Nplanes());
176 
177  fHitsByWire[0] = dir.make<TH1D>("HitsByWire0", ";Wire #", fGeometry->Nwires(0), 0., fGeometry->Nwires(0));
178  fHitsByWire[1] = dir.make<TH1D>("HitsByWire1", ";Wire #", fGeometry->Nwires(1), 0., fGeometry->Nwires(1));
179  fHitsByWire[2] = dir.make<TH1D>("HitsByWire2", ";Wire #", fGeometry->Nwires(2), 0., fGeometry->Nwires(2));
180 
181  fDriftTimes[0] = dir.make<TH1D>("DriftTime0", ";time(ticks)", 2048, 0., 4096.);
182  fDriftTimes[1] = dir.make<TH1D>("DriftTime1", ";time(ticks)", 2048, 0., 4096.);
183  fDriftTimes[2] = dir.make<TH1D>("DriftTime2", ";time(ticks)", 2048, 0., 4096.);
184 
185  fHitsByTime[0] = dir.make<TH1D>("HitsByTime0", ";Tick", 1024, 0., 4096.);
186  fHitsByTime[1] = dir.make<TH1D>("HitsByTime1", ";Tick", 1024, 0., 4096.);
187  fHitsByTime[2] = dir.make<TH1D>("HitsByTime2", ";Tick", 1024, 0., 4096.);
188 
189  fPulseHeight[0] = dir.make<TH1D>("PulseHeight0", "PH (ADC)", 300, 0., 150.);
190  fPulseHeight[1] = dir.make<TH1D>("PulseHeight1", "PH (ADC)", 300, 0., 150.);
191  fPulseHeight[2] = dir.make<TH1D>("PulseHeight2", "PH (ADC)", 300, 0., 150.);
192  fPulseHeightSingle[0] = dir.make<TH1D>("PulseHeightS0", "PH (ADC)", 300, 0., 150.);
193  fPulseHeightSingle[1] = dir.make<TH1D>("PulseHeightS1", "PH (ADC)", 300, 0., 150.);
194  fPulseHeightSingle[2] = dir.make<TH1D>("PulseHeightS2", "PH (ADC)", 300, 0., 150.);
195  fPulseHeightMultiMax[0] = dir.make<TH1D>("PulseHeightMM0", "PH (ADC)", 300, 0., 150.);
196  fPulseHeightMultiMax[1] = dir.make<TH1D>("PulseHeightMM1", "PH (ADC)", 300, 0., 150.);
197  fPulseHeightMultiMax[2] = dir.make<TH1D>("PulseHeightMM2", "PH (ADC)", 300, 0., 150.);
198  fPulseHeightMulti[0] = dir.make<TH1D>("PulseHeightM0", "PH (ADC)", 300, 0., 150.);
199  fPulseHeightMulti[1] = dir.make<TH1D>("PulseHeightM1", "PH (ADC)", 300, 0., 150.);
200  fPulseHeightMulti[2] = dir.make<TH1D>("PulseHeightM2", "PH (ADC)", 300, 0., 150.);
201  fPulseHeightLong[0] = dir.make<TH1D>("PulseHeightL0", "PH (ADC)", 300, 0., 150.);
202  fPulseHeightLong[1] = dir.make<TH1D>("PulseHeightL1", "PH (ADC)", 300, 0., 150.);
203  fPulseHeightLong[2] = dir.make<TH1D>("PulseHeightL2", "PH (ADC)", 300, 0., 150.);
204  fChi2DOF[0] = dir.make<TH1D>("Chi2DOF0", "Chi2DOF", 502, -1., 250.);
205  fChi2DOF[1] = dir.make<TH1D>("Chi2DOF1", "Chi2DOF", 502, -1., 250.);
206  fChi2DOF[2] = dir.make<TH1D>("Chi2DOF2", "Chi2DOF", 502, -1., 250.);
207  fNumDegFree[0] = dir.make<TH1D>("NumDegFree0", "NDF", 200, 0., 200.);
208  fNumDegFree[1] = dir.make<TH1D>("NumDegFree1", "NDF", 200, 0., 200.);
209  fNumDegFree[2] = dir.make<TH1D>("NumDegFree2", "NDF", 200, 0., 200.);
210  fChi2DOFSingle[0] = dir.make<TH1D>("Chi2DOFS0", "Chi2DOF", 502, -1., 250.);
211  fChi2DOFSingle[1] = dir.make<TH1D>("Chi2DOFS1", "Chi2DOF", 502, -1., 250.);
212  fChi2DOFSingle[2] = dir.make<TH1D>("Chi2DOFS2", "Chi2DOF", 502, -1., 250.);
213  fHitMult[0] = dir.make<TH1D>("HitMult0", "# hits", 15, 0., 15.);
214  fHitMult[1] = dir.make<TH1D>("HitMult1", "# hits", 15, 0., 15.);
215  fHitMult[2] = dir.make<TH1D>("HitMult2", "# hits", 15, 0., 15.);
216  fHitCharge[0] = dir.make<TH1D>("HitCharge0", "Charge", 1000, 0., 2000.);
217  fHitCharge[1] = dir.make<TH1D>("HitCharge1", "Charge", 1000, 0., 2000.);
218  fHitCharge[2] = dir.make<TH1D>("HitCharge2", "Charge", 1000, 0., 2000.);
219  fHitChargeSingle[0] = dir.make<TH1D>("HitChargeS0", "Charge", 1000, 0., 2000.);
220  fHitChargeSingle[1] = dir.make<TH1D>("HitChargeS1", "Charge", 1000, 0., 2000.);
221  fHitChargeSingle[2] = dir.make<TH1D>("HitChargeS2", "Charge", 1000, 0., 2000.);
222  fHitChargeMulti[0] = dir.make<TH1D>("HitChargeM0", "Charge", 1000, 0., 2000.);
223  fHitChargeMulti[1] = dir.make<TH1D>("HitChargeM1", "Charge", 1000, 0., 2000.);
224  fHitChargeMulti[2] = dir.make<TH1D>("HitChargeM2", "Charge", 1000, 0., 2000.);
225  fHitChargeLong[0] = dir.make<TH1D>("HitChargeL0", "Charge", 1000, 0., 2000.);
226  fHitChargeLong[1] = dir.make<TH1D>("HitChargeL1", "Charge", 1000, 0., 2000.);
227  fHitChargeLong[2] = dir.make<TH1D>("HitChargeL2", "Charge", 1000, 0., 2000.);
228  fFitWidth[0] = dir.make<TH1D>("FitWidth0", "Width", 100, 0., 20.);
229  fFitWidth[1] = dir.make<TH1D>("FitWidth1", "Width", 100, 0., 20.);
230  fFitWidth[2] = dir.make<TH1D>("FitWidth2", "Width", 100, 0., 20.);
231  fFitWidthSingle[0] = dir.make<TH1D>("FitWidthS0", "Width", 100, 0., 20.);
232  fFitWidthSingle[1] = dir.make<TH1D>("FitWidthS1", "Width", 100, 0., 20.);
233  fFitWidthSingle[2] = dir.make<TH1D>("FitWidthS2", "Width", 100, 0., 20.);
234  fFitWidthMulti[0] = dir.make<TH1D>("FitWidthM0", "Width", 100, 0., 20.);
235  fFitWidthMulti[1] = dir.make<TH1D>("FitWidthM1", "Width", 100, 0., 20.);
236  fFitWidthMulti[2] = dir.make<TH1D>("FitWidthM2", "Width", 100, 0., 20.);
237  fFitWidthLong[0] = dir.make<TH1D>("FitWidthL0", "Width", 100, 0., 20.);
238  fFitWidthLong[1] = dir.make<TH1D>("FitWidthL1", "Width", 100, 0., 20.);
239  fFitWidthLong[2] = dir.make<TH1D>("FitWidthL2", "Width", 100, 0., 20.);
240  fHitSumADC[0] = dir.make<TH1D>("SumADC0", "Sum ADC", 1000, 0., 2000.);
241  fHitSumADC[1] = dir.make<TH1D>("SumADC1", "Sum ADC", 1000, 0., 2000.);
242  fHitSumADC[2] = dir.make<TH1D>("SumADC2", "Sum ADC", 1000, 0., 2000.);
243 
244  fNDFVsChi2[0] = dir.make<TH2D>("NDFVsChi20", ";NDF;Chi2", 50, 0., 50., 101, -1., 100.);
245  fNDFVsChi2[1] = dir.make<TH2D>("NDFVsChi21", ";NDF;Chi2", 50, 0., 50., 101, -1., 100.);
246  fNDFVsChi2[2] = dir.make<TH2D>("NDFVsChi22", ";NDF;Chi2", 50, 0., 50., 101, -1., 100.);
247 
248  fPulseHVsWidth[0] = dir.make<TH2D>("PHVsWidth0", ";PH;Width", 100, 0., 100., 100, 0., 20.);
249  fPulseHVsWidth[1] = dir.make<TH2D>("PHVsWidth1", ";PH;Width", 100, 0., 100., 100, 0., 20.);
250  fPulseHVsWidth[2] = dir.make<TH2D>("PHVsWidth2", ";PH;Width", 100, 0., 100., 100, 0., 20.);
251 
252  fPulseHVsCharge[0] = dir.make<TH2D>("PHVsChrg0", ";PH;Q", 100, 0., 100., 100, 0., 2000.);
253  fPulseHVsCharge[1] = dir.make<TH2D>("PHVsChrg1", ";PH;Q", 100, 0., 100., 100, 0., 2000.);
254  fPulseHVsCharge[2] = dir.make<TH2D>("PHVsChrg2", ";PH;Q", 100, 0., 100., 100, 0., 2000.);
255 
256  fPulseHVsHitNo[0] = dir.make<TProfile>("PHVsNo0", ";Hit #;PH", 1000, 0., 1000., 0., 100.);
257  fPulseHVsHitNo[1] = dir.make<TProfile>("PHVsNo1", ";Hit #;PH", 1000, 0., 1000., 0., 100.);
258  fPulseHVsHitNo[2] = dir.make<TProfile>("PHVsNo2", ";Hit #;PH", 1000, 0., 1000., 0., 100.);
259 
260  fChargeVsHitNo[0] = dir.make<TProfile>("QVsNo0", ";Hit No;Q", 1000, 0., 1000., 0., 2000.);
261  fChargeVsHitNo[1] = dir.make<TProfile>("QVsNo1", ";Hit No;Q", 1000, 0., 1000., 0., 2000.);
262  fChargeVsHitNo[2] = dir.make<TProfile>("QVsNo2", ";Hit No;Q", 1000, 0., 1000., 0., 2000.);
263 
264  fChargeVsHitNoS[0] = dir.make<TProfile>("QVsNoS0", ";Hit No;Q", 1000, 0., 1000., 0., 2000.);
265  fChargeVsHitNoS[1] = dir.make<TProfile>("QVsNoS1", ";Hit No;Q", 1000, 0., 1000., 0., 2000.);
266  fChargeVsHitNoS[2] = dir.make<TProfile>("QVsNoS2", ";Hit No;Q", 1000, 0., 1000., 0., 2000.);
267 
268  fBadWPulseHeight = dir.make<TH1D>("BWPulseHeight", "PH (ADC)", 300, 0., 150.);
269  fBadWPulseHVsWidth = dir.make<TH2D>("BWPHVsWidth", ";PH;Width", 100, 0., 100., 100, 0., 10.);
270  fBadWHitsByWire = dir.make<TH1D>("BWHitsByWire", ";Wire #", fGeometry->Nwires(2), 0., fGeometry->Nwires(2));
271 
272  fSPHvsIdx[0] = dir.make<TH2D>("SPHVsIdx0", ";PH;Idx", 30, 0., 30., 100, 0., 100.);
273  fSPHvsIdx[1] = dir.make<TH2D>("SPHVsIdx1", ";PH;Idx", 30, 0., 30., 100, 0., 100.);
274  fSPHvsIdx[2] = dir.make<TH2D>("SPHVsIdx2", ";PH;Idx", 30, 0., 30., 100, 0., 100.);
275 
276  fSWidVsIdx[0] = dir.make<TH2D>("SWidsIdx0", ";Width;Idx", 30, 0., 30., 100, 0., 10.);
277  fSWidVsIdx[1] = dir.make<TH2D>("SWidsIdx1", ";Width;Idx", 30, 0., 30., 100, 0., 10.);
278  fSWidVsIdx[2] = dir.make<TH2D>("SWidsIdx2", ";Width;Idx", 30, 0., 30., 100, 0., 10.);
279 
280  f1PPHvsWid[0] = dir.make<TH2D>("1PPHVsWid0", ";PH;Width", 100, 0., 100., 100, 0., 10.);
281  f1PPHvsWid[1] = dir.make<TH2D>("1PPHVsWid1", ";PH;Width", 100, 0., 100., 100, 0., 10.);
282  f1PPHvsWid[2] = dir.make<TH2D>("1PPHVsWid2", ";PH;Width", 100, 0., 100., 100, 0., 10.);
283 
284  fSPPHvsWid[0] = dir.make<TH2D>("SPPHVsWid0", ";PH;Width", 100, 0., 100., 100, 0., 10.);
285  fSPPHvsWid[1] = dir.make<TH2D>("SPPHVsWid1", ";PH;Width", 100, 0., 100., 100, 0., 10.);
286  fSPPHvsWid[2] = dir.make<TH2D>("SPPHVsWid2", ";PH;Width", 100, 0., 100., 100, 0., 10.);
287 
288  fSOPHvsWid[0] = dir.make<TH2D>("SOPHVsWid0", ";PH;Width", 100, 0., 100., 100, 0., 10.);
289  fSOPHvsWid[1] = dir.make<TH2D>("SOPHVsWid1", ";PH;Width", 100, 0., 100., 100, 0., 10.);
290  fSOPHvsWid[2] = dir.make<TH2D>("SOPHVsWid2", ";PH;Width", 100, 0., 100., 100, 0., 10.);
291 
292  fPHRatVsIdx[0] = dir.make<TH2D>("PHRatVsIdx0", ";PHRat;Idx", 30, 0., 30., 51, 0., 1.02);
293  fPHRatVsIdx[1] = dir.make<TH2D>("PHRatVsIdx1", ";PHRat;Idx", 30, 0., 30., 51, 0., 1.02);
294  fPHRatVsIdx[2] = dir.make<TH2D>("PHRatVsIdx2", ";PHRat;Idx", 30, 0., 30., 51, 0., 1.02);
295 
296  return;
297 }
298 
299 void BasicHitAnalysis::fillHistograms(const HitPtrVec& hitPtrVec) const
300 {
301  // Keep track of number of hits per view
302  size_t nHitsPerView[] = {0,0,0};
303  size_t negCount(0);
304 
305  std::vector<const recob::Hit*> hitSnippetVec;
306 
307  // Loop the hits and make some plots
308  for(const auto& hitPtr : hitPtrVec)
309  {
310  // Extract interesting hit parameters
311  const geo::WireID& wireID = hitPtr->WireID();
312  float chi2DOF = std::min(hitPtr->GoodnessOfFit(),float(249.8));
313  int numDOF = hitPtr->DegreesOfFreedom();
314  int hitMult = hitPtr->Multiplicity();
315  float peakTime = std::min(float(3199.),hitPtr->PeakTime());
316  float charge = hitPtr->Integral();
317  float sumADC = hitPtr->SummedADC();
318  float hitPH = std::min(hitPtr->PeakAmplitude(),float(149.8));
319  float hitSigma = hitPtr->RMS();
320 
321  size_t plane = wireID.Plane;
322  size_t wire = wireID.Wire;
323 
324  if (hitPH < 0.)
325  {
326  negCount++;
327  //std::cout << "Hit Plane: " << plane << ", wire: " << wire << ", T: " << peakTime << ", PH: " << hitPH << ", charge: " << charge << ", sumADC: " << sumADC << std::endl;
328  }
329 
330  nHitsPerView[plane]++;
331 
332  fHitsByWire[plane]->Fill(wire,1.);
333  fHitsByTime[plane]->Fill(peakTime, 1.);
334  fPulseHeight[plane]->Fill(std::min(float(149.5),hitPH), 1.);
335  fChi2DOF[plane]->Fill(chi2DOF, 1.);
336  fNumDegFree[plane]->Fill(numDOF, 1.);
337  fHitMult[plane]->Fill(hitMult, 1.);
338  fHitCharge[plane]->Fill(charge, 1.);
339  fFitWidth[plane]->Fill(std::min(float(19.99),hitSigma), 1.);
340  fHitSumADC[plane]->Fill(sumADC, 1.);
341  fNDFVsChi2[plane]->Fill(numDOF, chi2DOF, 1.);
342  fDriftTimes[plane]->Fill(peakTime, 1.);
343 
344  if (hitMult == 1)
345  {
346  fPulseHeightSingle[plane]->Fill(std::min(float(149.5),hitPH), 1.);
347  fFitWidthSingle[plane]->Fill(std::min(float(19.99),hitSigma), 1.);
348  fHitChargeSingle[plane]->Fill(charge, 1.);
349  fChi2DOFSingle[plane]->Fill(chi2DOF, 1.);
350  fPulseHVsWidth[plane]->Fill(std::min(float(99.9),hitPH), std::min(float(19.99),hitSigma), 1.);
351  fPulseHVsCharge[plane]->Fill(std::min(float(99.9),hitPH), std::min(float(1999.),charge), 1.);
352 
353  if (plane == 2 && hitPH < 5 && hitSigma < 2.2)
354  {
355  std::cout << "++> Plane: " << plane << ", wire: " << wire << ", peakTime: " << peakTime << ", ph: " << hitPH << ", w: " << hitSigma << std::endl;
356 
357  fBadWPulseHeight->Fill(std::min(float(149.5),hitPH),1.);
358  fBadWPulseHVsWidth->Fill(std::min(float(99.9),hitPH), std::min(float(9.99),hitSigma), 1.);
359  fBadWHitsByWire->Fill(wire,1.);
360  }
361  }
362  else
363  {
364  if (numDOF > 1)
365  {
366  fPulseHeightMulti[plane]->Fill(std::min(float(149.5),hitPH), 1.);
367  fFitWidthMulti[plane]->Fill(std::min(float(19.99),hitSigma), 1.);
368  fHitChargeMulti[plane]->Fill(charge, 1.);
369  }
370  else
371  {
372  fPulseHeightLong[plane]->Fill(std::min(float(149.5),hitPH), 1.);
373  fFitWidthLong[plane]->Fill(std::min(float(19.99),hitSigma), 1.);
374  fHitChargeLong[plane]->Fill(charge, 1.);
375  }
376  }
377 
378  // Look at hits on snippets
379  if (!hitSnippetVec.empty() && hitSnippetVec.back()->LocalIndex() >= hitPtr->LocalIndex())
380  {
381  // Only worried about multi hit snippets
382  if (hitSnippetVec.size() > 1)
383  {
384  // Sort in order of largest to smallest pulse height
385  std::sort(hitSnippetVec.begin(),hitSnippetVec.end(),[](const auto* left, const auto* right){return left->PeakAmplitude() > right->PeakAmplitude();});
386 
387  float maxPulseHeight = hitSnippetVec.front()->PeakAmplitude();
388 
389  fPulseHeightMultiMax[plane]->Fill(std::min(float(149.5),maxPulseHeight), 1.);
390 
391  for(size_t idx = 0; idx < hitSnippetVec.size(); idx++)
392  {
393  float pulseHeight = hitSnippetVec.at(idx)->PeakAmplitude();
394  float pulseWid = hitSnippetVec.at(idx)->RMS();
395  float pulseHeightRatio = pulseHeight / maxPulseHeight;
396 
397  size_t plane = hitSnippetVec.at(idx)->WireID().Plane;
398 
399  fSPHvsIdx[plane]->Fill(idx, std::min(float(99.9),pulseHeight), 1.);
400  fSWidVsIdx[plane]->Fill(idx, std::min(float(9.99),pulseWid), 1.);
401  fPHRatVsIdx[plane]->Fill(idx, pulseHeightRatio, 1.);
402 
403  if (idx == 0) fSPPHvsWid[plane]->Fill(std::min(float(99.9),pulseHeight), std::min(float(9.99),pulseWid), 1.);
404  else fSOPHvsWid[plane]->Fill(std::min(float(99.9),pulseHeight), std::min(float(9.99),pulseWid), 1.);
405  }
406  }
407  else
408  {
409  float pulseHeight = hitSnippetVec.front()->PeakAmplitude();
410  float pulseWid = hitSnippetVec.front()->RMS();
411  size_t plane = hitSnippetVec.front()->WireID().Plane;
412 
413  f1PPHvsWid[plane]->Fill(std::min(float(99.9),pulseHeight), std::min(float(9.99),pulseWid), 1.);
414  }
415 
416  hitSnippetVec.clear();
417  }
418 
419  hitSnippetVec.push_back(hitPtr.get());
420  }
421 
422  return;
423 }
424 
425 // Useful for normalizing histograms
426 void BasicHitAnalysis::endJob(int numEvents)
427 {
428  if (!fHitsByWire.empty())
429  {
430  // Normalize wire profiles to be hits/event
431  double normFactor(1./numEvents);
432 
433  for(size_t idx = 0; idx < 3; idx++) fHitsByWire[idx]->Scale(normFactor);
434  }
435 
436  return;
437 }
438 
439 DEFINE_ART_CLASS_TOOL(BasicHitAnalysis)
440 }
BasicHitAnalysis::BasicHitAnalysis::configure
void configure(fhicl::ParameterSet const &pset) override
Definition: BasicHitAnalysis_tool.cc:163
BasicHitAnalysis::HitPtrVec
std::vector< art::Ptr< recob::Hit >> HitPtrVec
Definition: BasicHitAnalysis_tool.cc:44
ServiceUtil.h
Utilities related to art service access.
BasicHitAnalysis::BasicHitAnalysis::fillHistograms
void fillHistograms(const HitPtrVec &) const override
Interface for filling histograms.
Definition: BasicHitAnalysis_tool.cc:299
Hit.h
Declaration of signal hit object.
right
walls no right
Definition: selectors.fcl:105
BasicHitAnalysis::BasicHitAnalysis::fLocalDirName
std::string fLocalDirName
Fraction for truncated mean.
Definition: BasicHitAnalysis_tool.cc:89
BasicHitAnalysis::BasicHitAnalysis::BasicHitAnalysis
BasicHitAnalysis(fhicl::ParameterSet const &pset)
Constructor.
Definition: BasicHitAnalysis_tool.cc:141
BasicHitAnalysis::TrackViewHitMap
std::map< int, ViewHitMap > TrackViewHitMap
Definition: BasicHitAnalysis_tool.cc:46
geo::WireID::Wire
WireID_t Wire
Index of the wire within its plane.
Definition: geo_types.h:580
actions::Scale
Scale(size_t pos, T factor) -> Scale< T >
geo::GeometryCore::Nwires
unsigned int Nwires(unsigned int p, unsigned int tpc=0, unsigned int cstat=0) const
Returns the total number of wires in the specified plane.
Definition: GeometryCore.h:3345
geo::GeometryCore::Nplanes
unsigned int Nplanes(unsigned int tpc=0, unsigned int cstat=0) const
Returns the total number of wire planes in the specified TPC.
Definition: GeometryCore.h:2854
BasicHitAnalysis::BasicHitAnalysis::initializeHists
void initializeHists(art::ServiceHandle< art::TFileService > &, const std::string &) override
Interface for initializing the histograms to be filled.
Definition: BasicHitAnalysis_tool.cc:170
BasicHitAnalysis::BasicHitAnalysis::fGeometry
const geo::GeometryCore * fGeometry
pointer to Geometry service
Definition: BasicHitAnalysis_tool.cc:131
left
walls no left
Definition: selectors.fcl:105
geo::PlaneID::Plane
PlaneID_t Plane
Index of the plane within its TPC.
Definition: geo_types.h:493
geo::GeometryCore
Description of geometry of one entire detector.
Definition: GeometryCore.h:1476
dropbox.dir
tuple dir
Definition: dropbox.py:28
BasicHitAnalysis::ViewHitMap
std::map< size_t, HitPtrVec > ViewHitMap
Definition: BasicHitAnalysis_tool.cc:45
BasicHitAnalysis::BasicHitAnalysis::endJob
void endJob(int numEvents) override
Interface for method to executve at the end of run processing.
Definition: BasicHitAnalysis_tool.cc:426
geo::WireID::WireID
constexpr WireID()=default
Default constructor: an invalid TPC ID.
IHitHistogramTool::HitPtrVec
std::vector< art::Ptr< recob::Hit >> HitPtrVec
Interface for filling histograms.
Definition: IHitHistogramTool.h:55
tfs
art::ServiceHandle< art::TFileService > tfs
Definition: NuVertexChargeTree_module.cc:254
Geometry.h
art framework interface to geometry description
cout
BEGIN_PROLOG could also be cout
Definition: messageservice.fcl:10
Generated by   doxygen 1.8.5