All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
SpacePointAna_module.cc
Go to the documentation of this file.
1 //
2 // Name: SpacePointAna_module.cc
3 //
4 // Purpose: Module SpacePointAna.
5 //
6 //
7 // Configuration parameters.
8 //
9 // HitModuleLabel: Hit module label.
10 // UseClusterHits: If true, use clustered hits (otherwise all hits).
11 // ClusterModuleLabel: Cluster module label.
12 // UseMC: Use MC truth information.
13 // SpacePointAlgTime: SpacePointAlg configuration (loose time cuts).
14 // SpacePointAlgSep: SpacePointAlg configuration (loose separation cuts).
15 // SpacePointAlgDefault: SpacePointAlg configuration (default cuts).
16 //
17 // Created: 2-Aug-2011 H. Greenlee
18 //
19 
20 #include "art/Framework/Core/EDAnalyzer.h"
21 #include "art/Framework/Core/ModuleMacros.h"
22 #include "art/Framework/Principal/Event.h"
23 #include "art/Framework/Services/Registry/ServiceHandle.h"
24 #include "art_root_io/TFileService.h"
25 #include "canvas/Persistency/Common/FindManyP.h"
26 #include "messagefacility/MessageLogger/MessageLogger.h"
27 
28 #include "larcore/Geometry/Geometry.h"
29 #include "larcorealg/Geometry/PlaneGeo.h"
30 #include "larcorealg/Geometry/TPCGeo.h"
31 #include "lardata/ArtDataHelper/ToElement.h"
32 #include "lardata/DetectorInfoServices/DetectorClocksService.h"
33 #include "lardata/DetectorInfoServices/DetectorPropertiesService.h"
34 #include "lardataobj/RecoBase/Cluster.h"
35 #include "lardataobj/RecoBase/Hit.h"
36 #include "lardataobj/RecoBase/SpacePoint.h"
37 #include "larreco/RecoAlg/SpacePointAlg.h"
38 #include "larsim/MCCheater/BackTrackerService.h"
39 
40 #include "TH1F.h"
41 #include "TH2F.h"
42 
43 #include "range/v3/view.hpp"
44 
45 using lar::to_element;
46 using ranges::views::transform;
47 
48 namespace trkf {
49 
50  class SpacePointAna : public art::EDAnalyzer {
51  public:
52  explicit SpacePointAna(fhicl::ParameterSet const& pset);
53 
54  private:
55  void analyze(const art::Event& evt) override;
56 
57  void bookHistograms(bool mc);
58 
59  // Fcl Attributes.
60 
61  const SpacePointAlg fSptalgTime; // Algorithm object (increased time cut).
62  const SpacePointAlg fSptalgSep; // Algorithm object (increased sepataion cut).
63  const SpacePointAlg fSptalgDefault; // Algorithm object (default cuts).
64  std::string fHitModuleLabel;
65  bool fUseClusterHits;
66  std::string fClusterModuleLabel;
67  bool fUseMC;
68  double fMinX; // Minimum x.
69  double fMaxX; // Maximum x.
70  double fMinY; // Minimum y.
71  double fMaxY; // Maximum y.
72  double fMinZ; // Minimum z.
73  double fMaxZ; // Maximum z.
74 
75  // Histograms.
76 
77  bool fBooked; // Have histograms been booked yet?
78  TH1F* fHDTUE; // U-drift electrons time difference.
79  TH1F* fHDTVE; // V-drift electrons time difference.
80  TH1F* fHDTWE; // W-drift electrons time difference.
81  TH1F* fHDTUPull; // U-drift electrons time pull.
82  TH1F* fHDTVPull; // V-drift electrons time pull.
83  TH1F* fHDTWPull; // W-drift electrons time pull.
84  TH1F* fHDTUV; // U-V time difference.
85  TH1F* fHDTVW; // V-W time difference.
86  TH1F* fHDTWU; // W-U time difference.
87  TH2F* fHDTUVU; // U-V time difference vs. U.
88  TH2F* fHDTUVV; // U-V time difference vs. V.
89  TH2F* fHDTVWV; // V-W time difference vs. V.
90  TH2F* fHDTVWW; // V-W time difference vs. W.
91  TH2F* fHDTWUW; // W-U time difference vs. W.
92  TH2F* fHDTWUU; // W-U time difference vs. U.
93  TH1F* fHS; // Spatial separation.
94  TH1F* fHchisq; // Space point chisquare.
95  TH1F* fHx; // X position.
96  TH1F* fHy; // Y position.
97  TH1F* fHz; // Z position.
98  TH1F* fHAmpU; // U hit amplitude.
99  TH1F* fHAmpV; // V hit amplitude.
100  TH1F* fHAmpW; // W hit amplitude.
101  TH1F* fHAreaU; // U hit area.
102  TH1F* fHAreaV; // V hit area.
103  TH1F* fHAreaW; // W hit area.
104  TH1F* fHSumU; // U hit sum ADC.
105  TH1F* fHSumV; // V hit sum ADC.
106  TH1F* fHSumW; // W hit sum ADC.
107  TH1F* fHMCdx; // X residual (reco vs. mc truth).
108  TH1F* fHMCdy; // Y residual (reco vs. mc truth).
109  TH1F* fHMCdz; // Z residual (reco vs. mc truth).
110  TH1F* fHMCxpull; // X pull (reco vs. mc truth).
111  TH1F* fHMCypull; // Y pull (reco vs. mc truth).
112  TH1F* fHMCzpull; // Z pull (reco vs. mc truth).
113 
114  // Statistics.
115 
116  int fNumEvent;
117  };
118 
119  DEFINE_ART_MODULE(SpacePointAna)
120 
121  SpacePointAna::SpacePointAna(const fhicl::ParameterSet& pset)
122  //
123  // Purpose: Constructor.
124  //
125  // Arguments: pset - Module parameters.
126  //
127  : EDAnalyzer(pset)
128  , fSptalgTime(pset.get<fhicl::ParameterSet>("SpacePointAlgTime"))
129  , fSptalgSep(pset.get<fhicl::ParameterSet>("SpacePointAlgSep"))
130  , fSptalgDefault(pset.get<fhicl::ParameterSet>("SpacePointAlgDefault"))
131  , fHitModuleLabel(pset.get<std::string>("HitModuleLabel"))
132  , fUseClusterHits(pset.get<bool>("UseClusterHits"))
133  , fClusterModuleLabel(pset.get<std::string>("ClusterModuleLabel"))
134  , fUseMC(pset.get<bool>("UseMC"))
135  , fMinX(pset.get<double>("MinX", -1.e10))
136  , fMaxX(pset.get<double>("MaxX", 1.e10))
137  , fMinY(pset.get<double>("MinY", -1.e10))
138  , fMaxY(pset.get<double>("MaxY", 1.e10))
139  , fMinZ(pset.get<double>("MinZ", -1.e10))
140  , fMaxZ(pset.get<double>("MaxZ", 1.e10))
141  , fBooked(false)
142  , fHDTUE(0)
143  , fHDTVE(0)
144  , fHDTWE(0)
145  , fHDTUPull(0)
146  , fHDTVPull(0)
147  , fHDTWPull(0)
148  , fHDTUV(0)
149  , fHDTVW(0)
150  , fHDTWU(0)
151  , fHDTUVU(0)
152  , fHDTUVV(0)
153  , fHDTVWV(0)
154  , fHDTVWW(0)
155  , fHDTWUW(0)
156  , fHDTWUU(0)
157  , fHS(0)
158  , fHchisq(0)
159  , fHx(0)
160  , fHy(0)
161  , fHz(0)
162  , fHAmpU(0)
163  , fHAmpV(0)
164  , fHAmpW(0)
165  , fHAreaU(0)
166  , fHAreaV(0)
167  , fHAreaW(0)
168  , fHSumU(0)
169  , fHSumV(0)
170  , fHSumW(0)
171  , fHMCdx(0)
172  , fHMCdy(0)
173  , fHMCdz(0)
174  , fHMCxpull(0)
175  , fHMCypull(0)
176  , fHMCzpull(0)
177  , fNumEvent(0)
178  {
179  mf::LogInfo("SpacePointAna") << "SpacePointAna configured with the following parameters:\n"
180  << " HitModuleLabel = " << fHitModuleLabel << "\n"
181  << " UseClusterHits = " << fUseClusterHits << "\n"
182  << " ClusterModuleLabel = " << fClusterModuleLabel << "\n"
183  << " UseMC = " << fUseMC;
184  }
185 
186  void
187  SpacePointAna::bookHistograms(bool mc)
188  {
189  if (!fBooked) {
190  fBooked = true;
191 
192  art::ServiceHandle<geo::Geometry const> geom;
193  art::ServiceHandle<art::TFileService> tfs;
194  art::TFileDirectory dir = tfs->mkdir("sptana", "SpacePointAna histograms");
195 
196  unsigned int nwiresU = 0, nwiresV = 0, nwiresW = 0;
197 
198  // Figure out the number of wires in U, V, and W planes.
199 
200  // Loop over cryostats, tpcs, and planes.
201 
202  for (unsigned int cstat = 0; cstat < geom->Ncryostats(); ++cstat) {
203 
204  const geo::CryostatGeo& cryogeom = geom->Cryostat(cstat);
205  unsigned int const ntpc = cryogeom.NTPC();
206 
207  for (unsigned int tpc = 0; tpc < ntpc; ++tpc) {
208 
209  const geo::TPCGeo& tpcgeom = cryogeom.TPC(tpc);
210  unsigned int const nplane = tpcgeom.Nplanes();
211 
212  for (unsigned int plane = 0; plane < nplane; ++plane) {
213 
214  const geo::PlaneGeo& pgeom = tpcgeom.Plane(plane);
215  unsigned int nwires = pgeom.Nwires();
216  geo::View_t view = pgeom.View();
217  if (view == geo::kU)
218  nwiresU = nwires;
219  else if (view == geo::kV)
220  nwiresV = nwires;
221  else if (view == geo::kZ)
222  nwiresW = nwires;
223  }
224  }
225  }
226 
227  if (mc && fUseMC) {
228  fHDTUE = dir.make<TH1F>("MCDTUE", "U-Drift Electrons Time Difference", 100, -5., 5.);
229  fHDTVE = dir.make<TH1F>("MCDTVE", "V-Drift Electrons Time Difference", 100, -5., 5.);
230  fHDTWE = dir.make<TH1F>("MCDTWE", "W-Drift Electrons Time Difference", 100, -5., 5.);
231  fHDTUPull = dir.make<TH1F>("MCDTUPull", "U-Drift Electrons Time Pull", 100, -50., 50.);
232  fHDTVPull = dir.make<TH1F>("MCDTVPull", "V-Drift Electrons Time Pull", 100, -50., 50.);
233  fHDTWPull = dir.make<TH1F>("MCDTWPull", "W-Drift Electrons Time Pull", 100, -50., 50.);
234  }
235  if (!fSptalgTime.merge()) {
236  fHDTUV = dir.make<TH1F>("DTUV", "U-V time difference", 100, -20., 20.);
237  fHDTVW = dir.make<TH1F>("DTVW", "V-W time difference", 100, -20., 20.);
238  fHDTWU = dir.make<TH1F>("DTWU", "W-U time difference", 100, -20., 20.);
239  fHDTUVU = dir.make<TH2F>(
240  "DTUVU", "U-V time difference vs. U", 100, 0., double(nwiresU), 100, -20., 20.);
241  fHDTUVV = dir.make<TH2F>(
242  "DTUVV", "U-V time difference vs. V", 100, 0., double(nwiresV), 100, -20., 20.);
243  fHDTVWV = dir.make<TH2F>(
244  "DTVWV", "V-W time difference vs. V", 100, 0., double(nwiresV), 100, -20., 20.);
245  fHDTVWW = dir.make<TH2F>(
246  "DTVWW", "V-W time difference vs. W", 100, 0., double(nwiresW), 100, -20., 20.);
247  fHDTWUW = dir.make<TH2F>(
248  "DTWUW", "W-U time difference vs. W", 100, 0., double(nwiresW), 100, -20., 20.);
249  fHDTWUU = dir.make<TH2F>(
250  "DTWUU", "W-U time difference vs. U", 100, 0., double(nwiresU), 100, -20., 20.);
251  fHS = dir.make<TH1F>("DS", "Spatial Separation", 100, -2., 2.);
252  }
253  fHchisq = dir.make<TH1F>("chisq", "Chisquare", 100, 0., 20.);
254 
255  fHx = dir.make<TH1F>("xpos", "X Position", 100, 0., 2. * geom->DetHalfWidth());
256  fHy =
257  dir.make<TH1F>("ypos", "Y Position", 100, -geom->DetHalfHeight(), geom->DetHalfHeight());
258  fHz = dir.make<TH1F>("zpos", "Z Position", 100, 0., geom->DetLength());
259  fHAmpU = dir.make<TH1F>("ampU", "U Hit Amplitude", 50, 0., 50.);
260  fHAmpV = dir.make<TH1F>("ampV", "V Hit Amplitude", 50, 0., 50.);
261  fHAmpW = dir.make<TH1F>("ampW", "W Hit Amplitude", 50, 0., 50.);
262  fHAreaU = dir.make<TH1F>("areaU", "U Hit Area", 100, 0., 500.);
263  fHAreaV = dir.make<TH1F>("areaV", "V Hit Area", 100, 0., 500.);
264  fHAreaW = dir.make<TH1F>("areaW", "W Hit Area", 100, 0., 500.);
265  fHSumU = dir.make<TH1F>("sumU", "U Hit Sum ADC", 100, 0., 500.);
266  fHSumV = dir.make<TH1F>("sumV", "V Hit Sum ADC", 100, 0., 500.);
267  fHSumW = dir.make<TH1F>("sumW", "W Hit Sum ADC", 100, 0., 500.);
268  if (mc && fUseMC) {
269  fHMCdx = dir.make<TH1F>("MCdx", "X MC Residual", 100, -1., 1.);
270  fHMCdy = dir.make<TH1F>("MCdy", "Y MC Residual", 100, -1., 1.);
271  fHMCdz = dir.make<TH1F>("MCdz", "Z MC Residual", 100, -1., 1.);
272  fHMCxpull = dir.make<TH1F>("MCxpull", "X MC Pull", 100, -50., 50.);
273  fHMCypull = dir.make<TH1F>("MCypull", "Y MC Pull", 100, -50., 50.);
274  fHMCzpull = dir.make<TH1F>("MCzpull", "Z MC Pull", 100, -50., 50.);
275  }
276  }
277  }
278 
279  void
280  SpacePointAna::analyze(const art::Event& evt)
281  {
282  ++fNumEvent;
283 
284  // Make sure histograms are booked.
285 
286  bool mc = !evt.isRealData();
287  bookHistograms(mc);
288 
289  // Get Services.
290 
291  art::ServiceHandle<geo::Geometry const> geom;
292 
293  // Get Hits.
294 
295  art::PtrVector<recob::Hit> hits;
296 
297  if (fUseClusterHits) {
298 
299  // Get clusters.
300 
301  art::Handle<std::vector<recob::Cluster>> clusterh;
302  evt.getByLabel(fClusterModuleLabel, clusterh);
303 
304  // Get hits from all clusters.
305  art::FindManyP<recob::Hit> fm(clusterh, evt, fClusterModuleLabel);
306 
307  if (clusterh.isValid()) {
308  int nclus = clusterh->size();
309 
310  for (int i = 0; i < nclus; ++i) {
311  art::Ptr<recob::Cluster> pclus(clusterh, i);
312  std::vector<art::Ptr<recob::Hit>> clushits = fm.at(i);
313  int nhits = clushits.size();
314  hits.reserve(hits.size() + nhits);
315 
316  for (std::vector<art::Ptr<recob::Hit>>::const_iterator ihit = clushits.begin();
317  ihit != clushits.end();
318  ++ihit) {
319  hits.push_back(*ihit);
320  }
321  }
322  }
323  }
324  else {
325 
326  // Get unclustered hits.
327 
328  art::Handle<std::vector<recob::Hit>> hith;
329  evt.getByLabel(fHitModuleLabel, hith);
330  if (hith.isValid()) {
331  int nhits = hith->size();
332  hits.reserve(nhits);
333 
334  for (int i = 0; i < nhits; ++i)
335  hits.push_back(art::Ptr<recob::Hit>(hith, i));
336  }
337  }
338 
339  // Fill histograms that don't depend on space points.
340 
341  auto const clockData = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(evt);
342  auto const detProp =
343  art::ServiceHandle<detinfo::DetectorPropertiesService const>()->DataFor(evt, clockData);
344 
345  if (mc && fUseMC) {
346 
347  art::ServiceHandle<cheat::BackTrackerService const> bt_serv;
348 
349  // Loop over hits and fill hit-electron time difference histogram.
350 
351  for (auto const& hitPtr : hits) {
352  const recob::Hit& hit = *hitPtr;
353 
354  double tpeak = hit.PeakTime();
355  double terr = hit.SigmaPeakTime();
356 
357  bool tav_ok = true;
358  double tav = 0.;
359  try {
360  std::vector<double> hitxyz = bt_serv->HitToXYZ(clockData, hit);
361  tav = detProp.ConvertXToTicks(
362  hitxyz[0], hit.WireID().Plane, hit.WireID().TPC, hit.WireID().Cryostat);
363  }
364  catch (cet::exception& x) {
365  tav_ok = false;
366  }
367  if (tav_ok) {
368  if (hit.View() == geo::kU) {
369  fHDTUE->Fill(tpeak - tav);
370  fHDTUPull->Fill((tpeak - tav) / terr);
371  }
372  else if (hit.View() == geo::kV) {
373  fHDTVE->Fill(tpeak - tav);
374  fHDTVPull->Fill((tpeak - tav) / terr);
375  }
376  else if (hit.View() == geo::kZ) {
377  fHDTWE->Fill(tpeak - tav);
378  fHDTWPull->Fill((tpeak - tav) / terr);
379  }
380  else
381  throw cet::exception("SpacePointAna") << "Bad view = " << hit.View() << "\n";
382  }
383  }
384  }
385 
386  std::vector<recob::SpacePoint> spts1; // For time histograms.
387  std::vector<recob::SpacePoint> spts2; // For separation histogram.
388  std::vector<recob::SpacePoint> spts3; // Default cuts.
389 
390  // If nonzero time cut is specified, make space points using that
391  // time cut (for time histograms).
392 
393  if (!fSptalgTime.merge()) {
394  if (mc && fUseMC)
395  fSptalgTime.makeMCTruthSpacePoints(clockData, detProp, hits, spts1);
396  else
397  fSptalgTime.makeSpacePoints(clockData, detProp, hits, spts1);
398 
399  // Report number of space points.
400 
401  MF_LOG_DEBUG("SpacePointAna")
402  << "Found " << spts1.size() << " space points using special time cut.";
403  }
404 
405  // If nonzero separation cut is specified, make space points using that
406  // separation cut (for separation histogram).
407 
408  if (!fSptalgSep.merge()) {
409  if (mc && fUseMC)
410  fSptalgSep.makeMCTruthSpacePoints(clockData, detProp, hits, spts2);
411  else
412  fSptalgSep.makeSpacePoints(clockData, detProp, hits, spts2);
413 
414  // Report number of space points.
415 
416  MF_LOG_DEBUG("SpacePointAna")
417  << "Found " << spts2.size() << " space points using special seperation cut.";
418  }
419 
420  // Make space points using default cuts.
421 
422  if (mc && fUseMC)
423  fSptalgDefault.makeMCTruthSpacePoints(clockData, detProp, hits, spts3);
424  else
425  fSptalgDefault.makeSpacePoints(clockData, detProp, hits, spts3);
426 
427  // Report number of space points.
428 
429  MF_LOG_DEBUG("SpacePointAna") << "Found " << spts3.size()
430  << " space points using default cuts.";
431 
432  if (!fSptalgTime.merge()) {
433 
434  // Loop over space points and fill time histograms.
435 
436  for (auto const& spt : spts1) {
437  if (spt.XYZ()[0] < fMinX || spt.XYZ()[0] > fMaxX || spt.XYZ()[1] < fMinY ||
438  spt.XYZ()[1] > fMaxY || spt.XYZ()[2] < fMinZ || spt.XYZ()[2] > fMaxZ)
439  continue;
440 
441  // Get hits associated with this SpacePoint.
442 
443  const art::PtrVector<recob::Hit>& spthits = fSptalgTime.getAssociatedHits(spt);
444 
445  // Make a double loop over hits and fill hit time difference histograms.
446 
447  for (auto const& hit1 : spthits | transform(to_element)) {
448  geo::WireID hit1WireID = hit1.WireID();
449  unsigned int tpc1 = hit1WireID.TPC;
450  unsigned int plane1 = hit1WireID.Plane;
451  unsigned int wire1 = hit1WireID.Wire;
452 
453  geo::View_t view1 = hit1.View();
454  double t1 = fSptalgTime.correctedTime(detProp, hit1);
455 
456  for (auto const& hit2 : spthits | transform(to_element)) {
457  geo::WireID hit2WireID = hit2.WireID();
458  unsigned int tpc2 = hit2WireID.TPC;
459  unsigned int plane2 = hit2WireID.Plane;
460  unsigned int wire2 = hit2WireID.Wire;
461 
462  // Require same tpc, different view.
463 
464  if (tpc1 == tpc2 && plane1 != plane2) {
465 
466  geo::View_t view2 = hit2.View();
467  double t2 = fSptalgTime.correctedTime(detProp, hit2);
468 
469  if (view1 == geo::kU) {
470  if (view2 == geo::kV) {
471  fHDTUV->Fill(t1 - t2);
472  fHDTUVU->Fill(double(wire1), t1 - t2);
473  fHDTUVV->Fill(double(wire2), t1 - t2);
474  }
475  if (view2 == geo::kZ) {
476  fHDTWU->Fill(t2 - t1);
477  fHDTWUW->Fill(double(wire2), t2 - t1);
478  fHDTWUU->Fill(double(wire1), t2 - t1);
479  }
480  }
481  if (view1 == geo::kV) {
482  if (view2 == geo::kZ) {
483  fHDTVW->Fill(t1 - t2);
484  fHDTVWV->Fill(double(wire1), t1 - t2);
485  fHDTVWW->Fill(double(wire2), t1 - t2);
486  }
487  if (view2 == geo::kU) {
488  fHDTUV->Fill(t2 - t1);
489  fHDTUVU->Fill(double(wire2), t2 - t1);
490  fHDTUVV->Fill(double(wire1), t2 - t1);
491  }
492  }
493  if (view1 == geo::kZ) {
494  if (view2 == geo::kU) {
495  fHDTWU->Fill(t1 - t2);
496  fHDTWUW->Fill(double(wire1), t1 - t2);
497  fHDTWUU->Fill(double(wire2), t1 - t2);
498  }
499  if (view2 == geo::kV) {
500  fHDTVW->Fill(t2 - t1);
501  fHDTVWV->Fill(double(wire2), t2 - t1);
502  fHDTVWW->Fill(double(wire1), t2 - t1);
503  }
504  }
505  }
506  }
507  }
508  }
509 
510  // Loop over space points and fill seperation histograms.
511 
512  for (auto const& spt : spts2) {
513  if (spt.XYZ()[0] < fMinX || spt.XYZ()[0] > fMaxX || spt.XYZ()[1] < fMinY ||
514  spt.XYZ()[1] > fMaxY || spt.XYZ()[2] < fMinZ || spt.XYZ()[2] > fMaxZ)
515  continue;
516 
517  // Get hits associated with this SpacePoint.
518 
519  const art::PtrVector<recob::Hit>& spthits = fSptalgSep.getAssociatedHits(spt);
520 
521  // Fill separation histogram.
522 
523  double sep = fSptalgSep.separation(spthits);
524  fHS->Fill(sep);
525  }
526  }
527 
528  // Loop over default space points and fill histograms.
529 
530  for (auto const& spt : spts3) {
531  if (spt.XYZ()[0] < fMinX || spt.XYZ()[0] > fMaxX || spt.XYZ()[1] < fMinY ||
532  spt.XYZ()[1] > fMaxY || spt.XYZ()[2] < fMinZ || spt.XYZ()[2] > fMaxZ)
533  continue;
534 
535  fHchisq->Fill(spt.Chisq());
536  fHx->Fill(spt.XYZ()[0]);
537  fHy->Fill(spt.XYZ()[1]);
538  fHz->Fill(spt.XYZ()[2]);
539 
540  // Get hits associated with this SpacePoint.
541 
542  std::vector<art::Ptr<recob::Hit>> spthits;
543  const art::PtrVector<recob::Hit>& av_spthits = fSptalgDefault.getAssociatedHits(spt);
544  for (auto const& ptr : av_spthits) {
545  spthits.push_back(ptr);
546  }
547 
548  // Fill single hit histograms.
549 
550  for (auto const& hitPtr : spthits) {
551  const recob::Hit& hit = *hitPtr;
552 
553  geo::View_t view = hit.View();
554 
555  if (view == geo::kU) {
556  fHAmpU->Fill(hit.PeakAmplitude());
557  fHAreaU->Fill(hit.Integral());
558  fHSumU->Fill(hit.SummedADC());
559  }
560  if (view == geo::kV) {
561  fHAmpV->Fill(hit.PeakAmplitude());
562  fHAreaV->Fill(hit.Integral());
563  fHSumV->Fill(hit.SummedADC());
564  }
565  if (view == geo::kZ) {
566  fHAmpW->Fill(hit.PeakAmplitude());
567  fHAreaW->Fill(hit.Integral());
568  fHSumW->Fill(hit.SummedADC());
569  }
570  }
571 
572  if (mc && fUseMC) {
573 
574  art::ServiceHandle<cheat::BackTrackerService const> bt_serv;
575 
576  try {
577  std::vector<double> mcxyz = bt_serv->SpacePointHitsToWeightedXYZ(clockData, spthits);
578  fHMCdx->Fill(spt.XYZ()[0] - mcxyz[0]);
579  fHMCdy->Fill(spt.XYZ()[1] - mcxyz[1]);
580  fHMCdz->Fill(spt.XYZ()[2] - mcxyz[2]);
581  if (spt.ErrXYZ()[0] > 0.)
582  fHMCxpull->Fill((spt.XYZ()[0] - mcxyz[0]) / std::sqrt(spt.ErrXYZ()[0]));
583  if (spt.ErrXYZ()[2] > 0.)
584  fHMCypull->Fill((spt.XYZ()[1] - mcxyz[1]) / std::sqrt(spt.ErrXYZ()[2]));
585  if (spt.ErrXYZ()[5] > 0.)
586  fHMCzpull->Fill((spt.XYZ()[2] - mcxyz[2]) / std::sqrt(spt.ErrXYZ()[5]));
587  }
588  catch (cet::exception&) {
589  }
590  }
591  }
592  }
593 }
trkf::SpacePointAlg::merge
bool merge() const noexcept
Definition: SpacePointAlg.h:92
pmtana::std
double std(const std::vector< short > &wf, const double ped_mean, size_t start, size_t nsample)
Definition: UtilFunc.cxx:42
lar::to_element
constexpr to_element_t to_element
Definition: ToElement.h:24
icarus::waveform_operations::details::transform
static constexpr Sample_t transform(Sample_t sample)
Definition: WaveformOperations.h:145
x
process_name opflash particleana ie x
Definition: run_opflash_electron.fcl:90
geo::View_t
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
geo::kV
Planes which measure V.
Definition: geo_types.h:130
recob::Hit::WireID
geo::WireID WireID() const
Definition: Hit.h:233
geo::TPCGeo::Nplanes
unsigned int Nplanes() const
Number of planes in this tpc.
Definition: TPCGeo.h:165
Hit.h
Declaration of signal hit object.
get
get()
Definition: PMTconfigurationExtraction_module.cc:168
geo::TPCGeo
Geometry information for a single TPC.
Definition: TPCGeo.h:38
geo::CryostatID::Cryostat
CryostatID_t Cryostat
Index of cryostat.
Definition: geo_types.h:212
geo::kZ
Planes which measure Z direction.
Definition: geo_types.h:132
recob::Hit::Integral
float Integral() const
Integral under the calibrated signal waveform of the hit, in tick x ADC units.
Definition: Hit.h:224
recob::Hit::View
geo::View_t View() const
View for the plane of the hit.
Definition: Hit.h:232
geo::WireID::Wire
WireID_t Wire
Index of the wire within its plane.
Definition: geo_types.h:580
hit
process_name hit
Definition: cheaterreco.fcl:51
geo::CryostatGeo
Geometry information for a single cryostat.
Definition: CryostatGeo.h:43
trkf::SpacePointAlg::getAssociatedHits
const art::PtrVector< recob::Hit > & getAssociatedHits(const recob::SpacePoint &spt) const
Definition: SpacePointAlg.cxx:1495
trkf::SpacePointAna::fSptalgSep
const SpacePointAlg fSptalgSep
Definition: SpacePointAna_module.cc:62
recob::Hit::PeakAmplitude
float PeakAmplitude() const
The estimated amplitude of the hit at its peak, in ADC units.
Definition: Hit.h:221
trkf::SpacePointAlg::makeMCTruthSpacePoints
void makeMCTruthSpacePoints(detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, const art::PtrVector< recob::Hit > &hits, std::vector< recob::SpacePoint > &spts) const
Definition: SpacePointAlg.cxx:810
trkf::SpacePointAlg::separation
double separation(const art::PtrVector< recob::Hit > &hits) const
Definition: SpacePointAlg.cxx:191
lar::dump::vector
auto vector(Vector const &v)
Returns a manipulator which will print the specified array.
Definition: DumpUtils.h:265
geo::kU
Planes which measure U.
Definition: geo_types.h:129
geo::PlaneGeo::View
View_t View() const
Which coordinate does this plane measure.
Definition: PlaneGeo.h:184
trkf::SpacePointAlg::makeSpacePoints
void makeSpacePoints(detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, const art::PtrVector< recob::Hit > &hits, std::vector< recob::SpacePoint > &spts) const
Definition: SpacePointAlg.cxx:797
geo::PlaneGeo
Geometry information for a single wire plane.The plane is represented in the geometry by a solid whic...
Definition: PlaneGeo.h:82
trkf::SpacePointAna::analyze
void analyze(const art::Event &evt) override
Definition: SpacePointAna_module.cc:280
geo::CryostatGeo::NTPC
unsigned int NTPC() const
Number of TPCs in this cryostat.
Definition: CryostatGeo.h:181
geo::PlaneID::Plane
PlaneID_t Plane
Index of the plane within its TPC.
Definition: geo_types.h:493
trkf::SpacePointAlg::correctedTime
double correctedTime(detinfo::DetectorPropertiesData const &detProp, const recob::Hit &hit) const
Definition: SpacePointAlg.cxx:167
Cluster.h
Declaration of cluster object.
dropbox.dir
tuple dir
Definition: dropbox.py:28
recob::Hit::PeakTime
float PeakTime() const
Time of the signal peak, in tick units.
Definition: Hit.h:218
trkf::SpacePointAna::SpacePointAna
SpacePointAna(fhicl::ParameterSet const &pset)
Definition: SpacePointAna_module.cc:121
PlaneGeo.h
Encapsulate the construction of a single detector plane.
trkf::SpacePointAna::fSptalgTime
const SpacePointAlg fSptalgTime
Definition: SpacePointAna_module.cc:61
geo::CryostatGeo::TPC
const TPCGeo & TPC(unsigned int itpc) const
Return the itpc&#39;th TPC in the cryostat.
Definition: CryostatGeo.cxx:93
geo::PlaneGeo::Nwires
unsigned int Nwires() const
Number of wires in this plane.
Definition: PlaneGeo.h:269
recob::Hit::SummedADC
float SummedADC() const
The sum of calibrated ADC counts of the hit (0. by default)
Definition: Hit.h:223
trkf::SpacePointAna::fSptalgDefault
const SpacePointAlg fSptalgDefault
Definition: SpacePointAna_module.cc:63
geo::WireID::WireID
constexpr WireID()=default
Default constructor: an invalid TPC ID.
e10
then echo ***************************************echo Variable FHICL_FILE_PATH not found echo You porbably haven t set up larsoft echo Try setup uboonecode vXX_XX_XX q e10
Definition: find_fhicl.sh:6
recob::Hit::SigmaPeakTime
float SigmaPeakTime() const
Uncertainty for the signal peak, in tick units.
Definition: Hit.h:219
geo::TPCGeo::Plane
PlaneGeo const & Plane(geo::View_t view) const
Return the plane in the tpc with View_t view.
Definition: TPCGeo.cxx:263
recob::Hit
2D representation of charge deposited in the TDC/wire plane
Definition: Hit.h:48
tfs
art::ServiceHandle< art::TFileService > tfs
Definition: NuVertexChargeTree_module.cc:254
tca::evt
TCEvent evt
Definition: DataStructs.cxx:8
geo::TPCID::TPC
TPCID_t TPC
Index of the TPC within its cryostat.
Definition: geo_types.h:406
SpacePointAlg.h
Algorithm for generating space points from hits.
Geometry.h
art framework interface to geometry description
detsim::detProp
auto const detProp
Definition: SimWire_module.cc:157
TPCGeo.h
Encapsulate the construction of a single detector plane.
Generated by   doxygen 1.8.5