All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
SpacePointAnalysisMC_tool.cc
Go to the documentation of this file.
1 #include "icaruscode/Analysis/tools/IHitEfficiencyHistogramTool.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 #include "lardata/DetectorInfoServices/DetectorClocksService.h"
11 #include "canvas/Persistency/Common/FindManyP.h"
12 
13 #include "larcore/Geometry/Geometry.h"
14 #include "larcore/CoreUtils/ServiceUtil.h" // lar::providerFrom()
15 #include "lardata/DetectorInfoServices/DetectorClocksService.h"
16 #include "larevt/CalibrationDBI/Interface/ChannelStatusService.h"
17 #include "larevt/CalibrationDBI/Interface/ChannelStatusProvider.h"
18 
19 #include "lardataobj/RawData/RawDigit.h"
20 #include "lardataobj/RecoBase/Wire.h"
21 #include "lardataobj/RecoBase/Hit.h"
22 #include "lardataobj/RecoBase/SpacePoint.h"
23 #include "lardataobj/Simulation/SimChannel.h"
24 #include "lardataobj/Simulation/SimEnergyDeposit.h"
25 #include "nusimdata/SimulationBase/MCParticle.h"
26 
27 #include "larsim/Simulation/LArVoxelID.h"
28 
29 // Eigen
30 #include <Eigen/Dense>
31 
32 #include "TH1.h"
33 #include "TH2.h"
34 #include "TProfile.h"
35 #include "TProfile2D.h"
36 #include "TTree.h"
37 
38 #include <cmath>
39 #include <tuple>
40 #include <algorithm>
41 
42 namespace SpacePointAnalysisMC
43 {
44  ////////////////////////////////////////////////////////////////////////
45  //
46  // Class: SpacePointAnalysisMC
47  // Module Type: producer
48  // File: SpacePointAnalysisMC.cc
49  //
50  // The intent of this module is to provide methods for
51  // "analyzing" hits on waveforms
52  //
53  // Configuration parameters:
54  //
55  // TruncMeanFraction - the fraction of waveform bins to discard when
56  //
57  // Created by Tracy Usher (usher@slac.stanford.edu) on February 19, 2016
58  //
59  ////////////////////////////////////////////////////////////////////////
60 
61 // The following typedefs will, obviously, be useful
62 using HitPtrVec = std::vector<art::Ptr<recob::Hit>>;
63 using ViewHitMap = std::map<size_t,HitPtrVec>;
64 using TrackViewHitMap = std::map<int,ViewHitMap>;
65 
66 // Define object to keep track of hit/spacepoint related items
67 class HitSpacePointObj
68 {
69 public:
70  HitSpacePointObj() : fTree(nullptr) {}
71 
72  void setBranches(TTree* tree)
73  {
74  tree->Branch("NumIDEsHit0", "std::vector<int>", &fNumIDEsHit0Vec);
75  tree->Branch("NumIDEsHit1", "std::vector<int>", &fNumIDEsHit1Vec);
76  tree->Branch("NumIDEsHit2", "std::vector<int>", &fNumIDEsHit2Vec);
77  tree->Branch("NumIDEsSpacePoint", "std::vector<int>", &fNumIDEsSpacePointVec);
78 
79  tree->Branch("SPQuality", "std::vector<float>", &fSPQualityVec);
80  tree->Branch("SPTotalCharge", "std::vector<float>", &fSPTotalChargeVec);
81  tree->Branch("SPAsymmetry", "std::vector<float>", &fSPAsymmetryVec);
82  tree->Branch("SmallestPH", "std::vector<float>", &fSmallestPHVec);
83  tree->Branch("LargestPH", "std::vector<float>", &fLargestPHVec);
84  tree->Branch("AveragePH", "std::vector<float>", &fAveragePHVec);
85  tree->Branch("LargestDelT", "std::vector<float>", &fLargestDelTVec);
86 
87  tree->Branch("NumLongHitsSP", "std::vector<int>", &fNumLongHitsVec);
88  tree->Branch("NumPlanesSimMatch", "std::vector<int>", &fNumPlanesSimMatchVec);
89  tree->Branch("NumIntersectSet", "std::vector<int>", &fNumIntersectSetVec);
90 
91  tree->Branch("SimHitDeltaT0", "std::vector<int>", &fSimHitDeltaT0Vec);
92  tree->Branch("SimHitDeltaT1", "std::vector<int>", &fSimHitDeltaT1Vec);
93  tree->Branch("SimHitDeltaT2", "std::vector<int>", &fSimHitDeltaT2Vec);
94  tree->Branch("SimDeltaT10", "std::vector<int>", &fSimDelta10Vec);
95  tree->Branch("SimDeltaT11", "std::vector<int>", &fSimDelta21Vec);
96  tree->Branch("HitDeltaT10", "std::vector<int>", &fHitDelta10Vec);
97  tree->Branch("HitDeltaT11", "std::vector<int>", &fHitDelta21Vec);
98  tree->Branch("MaxElectronDep0", "std::vector<float>", &fBigElecDep0Vec);
99  tree->Branch("MaxElectronDep1", "std::vector<float>", &fBigElecDep1Vec);
100  tree->Branch("MaxElectronDep2", "std::vector<float>", &fBigElecDep2Vec);
101 
102  fTree = tree;
103  }
104 
105  void fill()
106  {
107  if (fTree) fTree->Fill();
108  }
109 
110  void clear()
111  {
112  fNumIDEsHit0Vec.clear();
113  fNumIDEsHit1Vec.clear();
114  fNumIDEsHit2Vec.clear();
115  fNumIDEsSpacePointVec.clear();
116 
117  fSPQualityVec.clear();
118  fSPTotalChargeVec.clear();
119  fSPAsymmetryVec.clear();
120  fSmallestPHVec.clear();
121  fLargestPHVec.clear();
122  fAveragePHVec.clear();
123  fLargestDelTVec.clear();
124 
125  fNumLongHitsVec.clear();
126  fNumPlanesSimMatchVec.clear();
127  fNumIntersectSetVec.clear();
128 
129  fSimHitDeltaT0Vec.clear();
130  fSimHitDeltaT1Vec.clear();
131  fSimHitDeltaT2Vec.clear();
132  fSimDelta10Vec.clear();
133  fSimDelta21Vec.clear();
134  fHitDelta10Vec.clear();
135  fHitDelta21Vec.clear();
136  fBigElecDep0Vec.clear();
137  fBigElecDep1Vec.clear();
138  fBigElecDep2Vec.clear();
139  }
140 
141  // Define tuple vars, make public for direct access
142  std::vector<int> fNumIDEsHit0Vec;
143  std::vector<int> fNumIDEsHit1Vec;
144  std::vector<int> fNumIDEsHit2Vec;
145  std::vector<int> fNumIDEsSpacePointVec;
146 
147  std::vector<float> fSPQualityVec;
148  std::vector<float> fSPTotalChargeVec;
149  std::vector<float> fSPAsymmetryVec;
150  std::vector<float> fSmallestPHVec;
151  std::vector<float> fLargestPHVec;
152  std::vector<float> fAveragePHVec;
153  std::vector<float> fLargestDelTVec;
154 
155  std::vector<int> fNumLongHitsVec;
156  std::vector<int> fNumPlanesSimMatchVec;
157  std::vector<int> fNumIntersectSetVec;
158 
159  std::vector<int> fSimHitDeltaT0Vec;
160  std::vector<int> fSimHitDeltaT1Vec;
161  std::vector<int> fSimHitDeltaT2Vec;
162  std::vector<int> fSimDelta10Vec;
163  std::vector<int> fSimDelta21Vec;
164  std::vector<int> fHitDelta10Vec;
165  std::vector<int> fHitDelta21Vec;
166  std::vector<float> fBigElecDep0Vec;
167  std::vector<float> fBigElecDep1Vec;
168  std::vector<float> fBigElecDep2Vec;
169 
170 private:
171  TTree* fTree;
172 };
173 
174 
175 // Define object to keep track of hit/sim related tuple items
176 class HitSimulationTupleObj
177 {
178 public:
179  HitSimulationTupleObj() : fTree(nullptr) {}
180 
181  void setBranches(TTree* tree)
182  {
183  tree->Branch("TicksSimChannel", "std::vector<int>", &fTicksSimChannelVec);
184  tree->Branch("TicksSimChanMost", "std::vector<int>", &fTicksSimChanMostVec);
185  tree->Branch("TicksTotHit", "std::vector<int>", &fTicksTotHitVec);
186  tree->Branch("TicksMaxSimRel", "std::vector<int>", &fTicksMaxSimRelVec);
187  tree->Branch("TicksDiffSimHit", "std::vector<int>", &fTicksDiffSimHitVec);
188  tree->Branch("EneTotDepHit", "std::vector<float>", &fEneTotDepHitVec);
189  tree->Branch("NElecTotalHit", "std::vector<float>", &fNElecTotHitVec); //< Total number elecrons (all sources) for hit
190  tree->Branch("EneBestDepHit", "std::vector<float>", &fEneBestDepHitVec);
191  tree->Branch("NElecBestHit", "std::vector<float>", &fNElecBestHitVec); //< # electrons from primary track
192  tree->Branch("EneMaxDepHit", "std::vector<float>", &fEneMaxDepHitVec);
193  tree->Branch("NDF", "std::vector<int>", &fNDFHitVec); //< Number of degrees of freedom of hit fit
194  tree->Branch("Multiplicity", "std::vector<int>", &fMultiplicityHitVec); //< Multiplicity of the snippet the hit is on
195  tree->Branch("LocalIndex", "std::vector<int>", &fLocalIndexHitVec); //< The index of the hit within the snippet
196  tree->Branch("TimeOrder", "std::vector<int>", &fTimeOrderHitVec); //< Time order of the hit (selection variable)
197  tree->Branch("ChiSquare", "std::vector<float>", &fChiSquareHitVec); //< Chi square of fit
198  tree->Branch("SummedADC", "std::vector<float>", &fSummedADCHitVec); //< Sum of all ADC values start/end of snippet
199  tree->Branch("Integral", "std::vector<float>", &fIntegralHitVec); //< Integrated charge +/- n sigma about peak center
200  tree->Branch("PulseHeight", "std::vector<float>", &fPHHitVec); //< Pulse height of hit
201  tree->Branch("RMS", "std::vector<float>", &fRMSHitVec); //< RMS of hit (from fit)
202  tree->Branch("PHFraction", "std::vector<float>", &fPHFractionVec); //< Fraction of max pulse height this snippet
203 
204  tree->Branch("PulseHeightOrder", "std::vector<int>", &fPHOrderHitVec); //< Local index ordered by pulse height
205 
206  fTree = tree;
207 
208  return;
209  }
210 
211  void fill()
212  {
213  if (fTree) fTree->Fill();
214  }
215 
216  void clear()
217  {
218  fTicksSimChannelVec.clear();
219  fTicksSimChanMostVec.clear();
220  fTicksTotHitVec.clear();
221  fTicksMaxSimRelVec.clear();
222  fTicksDiffSimHitVec.clear();
223  fEneTotDepHitVec.clear();
224  fNElecTotHitVec.clear();
225  fEneBestDepHitVec.clear();
226  fNElecBestHitVec.clear();
227  fEneMaxDepHitVec.clear();
228  fNDFHitVec.clear();
229  fMultiplicityHitVec.clear();
230  fLocalIndexHitVec.clear();
231  fTimeOrderHitVec.clear();
232  fChiSquareHitVec.clear();
233  fSummedADCHitVec.clear();
234  fIntegralHitVec.clear();
235  fPHHitVec.clear();
236  fRMSHitVec.clear();
237  fPHFractionVec.clear();
238 
239  fPHOrderHitVec.clear();
240  }
241 
242  void fillSimInfo(int ticksSimChannel,
243  int ticksSimChanMost,
244  float totDepEne,
245  float totNumElectrons,
246  float bestDepEne,
247  float bestNumElectrons,
248  float maxDepEneTick
249  )
250  {
251  fTicksSimChannelVec.emplace_back(ticksSimChannel);
252  fTicksSimChanMostVec.emplace_back(ticksSimChanMost);
253  fEneTotDepHitVec.emplace_back(totDepEne);
254  fNElecTotHitVec.emplace_back(totNumElectrons);
255  fEneBestDepHitVec.emplace_back(bestDepEne);
256  fNElecBestHitVec.emplace_back(bestNumElectrons);
257  fEneMaxDepHitVec.emplace_back(maxDepEneTick);
258  }
259 
260  void fillMixedInfo(int hitWidth,
261  int ticksToMax,
262  int deltaTicks)
263  {
264  fTicksTotHitVec.emplace_back(hitWidth);
265  fTicksMaxSimRelVec.emplace_back(ticksToMax);
266  fTicksDiffSimHitVec.emplace_back(deltaTicks);
267  }
268 
269  void fillHitInfo(const recob::Hit* hit,int hitOrder)
270  {
271  fNDFHitVec.emplace_back(hit->DegreesOfFreedom());
272  fMultiplicityHitVec.emplace_back(hit->Multiplicity());
273  fLocalIndexHitVec.emplace_back(hit->LocalIndex());
274  fTimeOrderHitVec.emplace_back(hitOrder);
275  fChiSquareHitVec.emplace_back(hit->GoodnessOfFit());
276  fSummedADCHitVec.emplace_back(hit->SummedADC());
277  fIntegralHitVec.emplace_back(hit->Integral());
278  fPHHitVec.emplace_back(hit->PeakAmplitude());
279  fRMSHitVec.emplace_back(hit->RMS());
280  }
281 
282  // Define tuple values, these are public so can be diretly accessed for filling
283  std::vector<int> fTicksSimChannelVec;
284  std::vector<int> fTicksSimChanMostVec;
285  std::vector<int> fTicksTotHitVec;
286  std::vector<int> fTicksMaxSimRelVec;
287  std::vector<int> fTicksDiffSimHitVec;
288  std::vector<float> fEneTotDepHitVec;
289  std::vector<float> fNElecTotHitVec;
290  std::vector<float> fEneBestDepHitVec;
291  std::vector<float> fNElecBestHitVec;
292  std::vector<float> fEneMaxDepHitVec;
293  std::vector<int> fNDFHitVec;
294  std::vector<int> fMultiplicityHitVec;
295  std::vector<int> fLocalIndexHitVec;
296  std::vector<int> fTimeOrderHitVec;
297  std::vector<float> fChiSquareHitVec;
298  std::vector<float> fSummedADCHitVec;
299  std::vector<float> fIntegralHitVec;
300  std::vector<float> fPHHitVec;
301  std::vector<float> fRMSHitVec;
302 
303  std::vector<float> fPHFractionVec;
304  std::vector<int> fPHOrderHitVec;
305 
306 private:
307  TTree* fTree;
308 };
309 
310 
311 class SpacePointAnalysisMC : virtual public IHitEfficiencyHistogramTool
312 {
313 public:
314  /**
315  * @brief Constructor
316  *
317  * @param pset
318  */
319  explicit SpacePointAnalysisMC(fhicl::ParameterSet const & pset);
320 
321  /**
322  * @brief Destructor
323  */
324  ~SpacePointAnalysisMC();
325 
326  // provide for initialization
327  void configure(fhicl::ParameterSet const & pset) override;
328 
329  /**
330  * @brief Interface for initializing the histograms to be filled
331  *
332  * @param TFileService handle to the TFile service
333  * @param string subdirectory to store the hists in
334  */
335  void initializeHists(art::ServiceHandle<art::TFileService>&, const std::string&) override;
336 
337  /**
338  * @brief Interface for initializing the tuple variables
339  *
340  * @param TTree pointer to a TTree object to which to add variables
341  */
342  void initializeTuple(TTree*) override;
343 
344  /**
345  * @brief Interface for method to executve at the end of run processing
346  *
347  * @param int number of events to use for normalization
348  */
349  void endJob(int numEvents) override;
350 
351  /**
352  * @brief Interface for filling histograms
353  */
354  void fillHistograms(const art::Event&) const override;
355 
356 private:
357 
358  // Clear mutable variables
359  void clear() const;
360 
361  // Create a struct allowing us to sort IDEs in a set by largest to smallest energy
362  struct ideCompare
363  {
364  bool operator() (const sim::IDE* left, const sim::IDE* right) const {return left->energy > right->energy;}
365  };
366 
367  // Define structures for relating SimChannel to Voxels
368  using SimIDESet = std::set<const sim::IDE*,ideCompare>;
369  using IDEToVoxelIDMap = std::unordered_map<const sim::IDE*, sim::LArVoxelID>;
370  using VoxelIDToIDESetMap = std::map<sim::LArVoxelID, SimIDESet>;
371  using TDCToIDEMap = std::map<unsigned short, SimIDESet>; // We need this one in order
372  using ChanToTDCToIDEMap = std::map<raw::ChannelID_t, TDCToIDEMap>;
373  using VoxelIDSet = std::set<sim::LArVoxelID>;
374 
375  // And, of course, what we need is to be able to track a voxel back to the IDEs in each tick on each plane
376  using TDCToIDESetMap = std::unordered_map<unsigned short,SimIDESet>;
377  using PlaneToTDCToIDESetMap = std::map<unsigned short, TDCToIDESetMap>;
378  using VoxelIDToPlaneTDCIDEMap = std::map<sim::LArVoxelID, PlaneToTDCToIDESetMap>;
379 
380  // The following creates a trackID mapping
381  using TDCIDEPair = std::pair<unsigned short, const sim::IDE*>;
382  using TickTDCIDEVec = std::vector<TDCIDEPair>;
383  using ChanToTDCIDEMap = std::unordered_map<raw::ChannelID_t,TickTDCIDEVec>;
384  using TrackIDChanToTDCIDEMap = std::unordered_map<int,ChanToTDCIDEMap>;
385 
386  // More data structures, here we want to keep track of the start/peak/end of the charge deposit along a wire for a given track
387  using ChargeDeposit = std::tuple<TDCIDEPair,TDCIDEPair,TDCIDEPair,float,float>;
388  using ChargeDepositVec = std::vector<ChargeDeposit>;
389  using ChanToChargeMap = std::map<raw::ChannelID_t,ChargeDepositVec>;
390  using TrackToChanChargeMap = std::unordered_map<int,ChanToChargeMap>;
391 
392  // Define a function to map IDE's from SimChannel objects to Track IDs
393  void makeTrackToChanChargeMap(const TrackIDChanToTDCIDEMap&, TrackToChanChargeMap&, float&, int&) const;
394 
395  // Relate hits to voxels
396  using HitPointerVec = std::vector<const recob::Hit*>;
397  using RecobHitToVoxelIDMap = std::unordered_map<const recob::Hit*, VoxelIDSet>;
398 
399  void compareHitsToSim(const art::Event&, const ChanToTDCToIDEMap&, const ChanToChargeMap&, const ChanToTDCIDEMap&, const IDEToVoxelIDMap&, RecobHitToVoxelIDMap&) const;
400 
401  void matchHitSim(const detinfo::DetectorClocksData& clockData,
402  const HitPointerVec&, const ChanToTDCToIDEMap&, const ChargeDepositVec&, const ChanToTDCIDEMap&, const IDEToVoxelIDMap&, RecobHitToVoxelIDMap&) const;
403 
404  void compareSpacePointsToSim(const art::Event&,
405  const detinfo::DetectorClocksData& clockData,
406  const RecobHitToVoxelIDMap&, const VoxelIDToPlaneTDCIDEMap&) const;
407 
408  // Fcl parameters.
409  std::vector<art::InputTag> fRecobHitLabelVec;
410  std::vector<art::InputTag> fSpacePointLabelVec;
411  art::InputTag fMCParticleProducerLabel;
412  art::InputTag fSimChannelProducerLabel;
413  art::InputTag fSimEnergyProducerLabel;
414  art::InputTag fBadChannelProducerLabel;
415  std::vector<int> fOffsetVec; ///< Allow offsets for each plane
416  float fSimChannelMinEnergy;
417  float fSimEnergyMinEnergy;
418 
419  // TTree variables
420  mutable TTree* fTree;
421 
422  mutable std::vector<int> fTPCVec;
423  mutable std::vector<int> fCryoVec;
424  mutable std::vector<int> fPlaneVec;
425 
426  using HitSimObjVec = std::vector<HitSimulationTupleObj>;
427 
428  mutable HitSimObjVec fHitSimObjVec;
429  mutable HitSpacePointObj fHitSpacePointObj;
430 
431  // Useful services, keep copies for now (we can update during begin run periods)
432  const geo::GeometryCore* fGeometry; ///< pointer to Geometry service
433 };
434 
435 //----------------------------------------------------------------------------
436 /// Constructor.
437 ///
438 /// Arguments:
439 ///
440 /// pset - Fcl parameters.
441 ///
442 SpacePointAnalysisMC::SpacePointAnalysisMC(fhicl::ParameterSet const & pset) : fTree(nullptr)
443 {
444  fGeometry = lar::providerFrom<geo::Geometry>();
445 
446  configure(pset);
447 
448  // Report.
449  mf::LogInfo("SpacePointAnalysisMC") << "SpacePointAnalysisMC configured\n";
450 }
451 
452 //----------------------------------------------------------------------------
453 /// Destructor.
454 SpacePointAnalysisMC::~SpacePointAnalysisMC()
455 {}
456 
457 //----------------------------------------------------------------------------
458 /// Reconfigure method.
459 ///
460 /// Arguments:
461 ///
462 /// pset - Fcl parameter set.
463 ///
464 void SpacePointAnalysisMC::configure(fhicl::ParameterSet const & pset)
465 {
466  fRecobHitLabelVec = pset.get< std::vector<art::InputTag>>("HitLabelVec", std::vector<art::InputTag>() = {"cluster3d"});
467  fSpacePointLabelVec = pset.get< std::vector<art::InputTag>>("SpacePointLabelVec", std::vector<art::InputTag>() = {"cluster3d"});
468  fMCParticleProducerLabel = pset.get< art::InputTag >("MCParticleLabel", "largeant");
469  fSimChannelProducerLabel = pset.get< art::InputTag >("SimChannelLabel", "largeant");
470  fSimEnergyProducerLabel = pset.get< art::InputTag >("SimEnergyLabel", "largeant");
471  fOffsetVec = pset.get<std::vector<int> >("OffsetVec", std::vector<int>()={0,0,0});
472  fSimChannelMinEnergy = pset.get<float >("SimChannelMinEnergy", std::numeric_limits<float>::epsilon());
473  fSimEnergyMinEnergy = pset.get<float >("SimEnergyMinEnergy", std::numeric_limits<float>::epsilon());
474 
475  return;
476 }
477 
478 //----------------------------------------------------------------------------
479 /// Begin job method.
480 void SpacePointAnalysisMC::initializeHists(art::ServiceHandle<art::TFileService>& tfs, const std::string& dirName)
481 {
482  // Make a directory for these histograms
483 // art::TFileDirectory dir = tfs->mkdir(dirName.c_str());
484 
485  return;
486 }
487 
488 void SpacePointAnalysisMC::initializeTuple(TTree* tree)
489 {
490  // Access ART's TFileService, which will handle creating and writing
491  // histograms and n-tuples for us.
492  art::ServiceHandle<art::TFileService> tfs;
493 
494  fTree = tree;
495 
496  fTree->Branch("CryostataVec", "std::vector<int>", &fCryoVec);
497  fTree->Branch("TPCVec", "std::vector<int>", &fTPCVec);
498  fTree->Branch("PlaneVec", "std::vector<int>", &fPlaneVec);
499 
500  // Set up specific branch for space points
501  TTree* locTree = tfs->makeAndRegister<TTree>("SpacePoint_t","SpacePoint Tuple");
502 
503  fHitSpacePointObj.setBranches(locTree);
504 
505  fHitSimObjVec.resize(fGeometry->Nplanes());
506 
507  for(size_t plane = 0; plane < fGeometry->Nplanes(); plane++)
508  {
509  // Set up specific branch for space points
510  locTree = tfs->makeAndRegister<TTree>("MatchedHits_P"+std::to_string(plane),"Matched Hits Tuple plane "+std::to_string(plane));
511 
512  fHitSimObjVec[plane].setBranches(locTree);
513  }
514 
515  clear();
516 
517  return;
518 }
519 
520 void SpacePointAnalysisMC::clear() const
521 {
522  fTPCVec.clear();
523  fCryoVec.clear();
524  fPlaneVec.clear();
525 
526  fHitSpacePointObj.clear();
527 
528  for(auto& hitObj : fHitSimObjVec) hitObj.clear();
529 
530  return;
531 }
532 
533 void SpacePointAnalysisMC::fillHistograms(const art::Event& event) const
534 {
535  // Ok... this is starting to grow too much and get out of control... we will need to break it up directly...
536 
537  // Always clear the tuple
538  clear();
539 
540  art::Handle<std::vector<sim::SimChannel>> simChannelHandle;
541  event.getByLabel(fSimChannelProducerLabel, simChannelHandle);
542 
543  if (!simChannelHandle.isValid() || simChannelHandle->empty() ) return;
544 
545  art::Handle<std::vector<sim::SimEnergyDeposit>> simEnergyHandle;
546  event.getByLabel(fSimEnergyProducerLabel, simEnergyHandle);
547 
548 // if (!simEnergyHandle.isValid() || simEnergyHandle->empty()) return;
549 
550  art::Handle<std::vector<simb::MCParticle>> mcParticleHandle;
551  event.getByLabel(fMCParticleProducerLabel, mcParticleHandle);
552 
553  // If there is no sim channel informaton then exit
554  if (!mcParticleHandle.isValid()) return;
555 
556  mf::LogDebug("SpacePointAnalysisMC") << "+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++" << std::endl;
557 
558  // First task is to build a map between ides and voxel ids (that we calcualate based on position)
559  // and also get the reverse since it will be useful in the end.
560  // At the same time should also build a mapping of ides per channel so we can do quick hit lookup
561  IDEToVoxelIDMap ideToVoxelIDMap;
562  VoxelIDToIDESetMap voxelIDToIDEMap;
563  ChanToTDCToIDEMap chanToTDCToIDEMap;
564  VoxelIDSet simChannelVoxelIDSet;
565  VoxelIDToPlaneTDCIDEMap voxelIDToPlaneTDCIDEMap;
566 
567  TrackIDChanToTDCIDEMap trackIDChanToTDCIDEMap;
568 
569  // Fill the above maps/structures
570  for(const auto& simChannel : *simChannelHandle)
571  {
572  raw::ChannelID_t channel = simChannel.Channel();
573 
574  geo::WireID wireID = fGeometry->ChannelToWire(channel).front();
575 
576  for(const auto& tdcide : simChannel.TDCIDEMap())
577  {
578  for(const auto& ide : tdcide.second) //chanToTDCToIDEMap[simChannel.Channel()][tdcide.first] = ide;
579  {
580  if (ide.energy < fSimChannelMinEnergy) continue;
581 
582  sim::LArVoxelID voxelID(ide.x,ide.y,ide.z,0.);
583 
584  ideToVoxelIDMap[&ide] = voxelID;
585  voxelIDToIDEMap[voxelID].insert(&ide);
586  chanToTDCToIDEMap[simChannel.Channel()][tdcide.first].insert(&ide);
587  simChannelVoxelIDSet.insert(voxelID);
588 
589  trackIDChanToTDCIDEMap[ide.trackID][simChannel.Channel()].emplace_back(tdcide.first,&ide);
590 
591  voxelIDToPlaneTDCIDEMap[voxelID][wireID.Plane][tdcide.first].insert(&ide);
592 
593  if (ide.energy < std::numeric_limits<float>::epsilon()) mf::LogDebug("SpacePointAnalysisMC") << ">> epsilon simchan deposited energy: " << ide.energy << std::endl;
594  }
595  }
596  }
597 
598  // More data structures, here we want to keep track of the start/peak/end of the charge deposit along a wire for a given track
599  TrackToChanChargeMap trackToChanChargeMap;
600 
601  // Go through the list of track to ides and get the total deposited energy per track
602  float bestTotDepEne(0.);
603  int bestTrackID(0);
604 
605  makeTrackToChanChargeMap(trackIDChanToTDCIDEMap, trackToChanChargeMap, bestTotDepEne, bestTrackID);
606 
607  // Ok, for my next trick I want to build a mapping between hits and voxel IDs. Note that any given hit can be associated to more than one voxel...
608  // We do this on the entire hit collection, ultimately we will want to consider SpacePoint efficiency (this could be done in the loop over SpacePoints
609  // using the associated hits and would save time/memory)
610  using VoxelIDSet = std::set<sim::LArVoxelID>;
611  using RecobHitToVoxelIDMap = std::unordered_map<const recob::Hit*, VoxelIDSet>;
612 
613  RecobHitToVoxelIDMap recobHitToVoxelIDMap;
614 
615  ChanToTDCIDEMap& chanToTDCIDEMap = trackIDChanToTDCIDEMap[bestTrackID];
616 
617  // Recover the "best" track info to start
618  TrackToChanChargeMap::const_iterator chanToChargeMapItr = trackToChanChargeMap.find(bestTrackID);
619 
620  // Process the hit/simulation
621  compareHitsToSim(event, chanToTDCToIDEMap, chanToChargeMapItr->second, chanToTDCIDEMap, ideToVoxelIDMap, recobHitToVoxelIDMap);
622 
623  // Now do the space points
624  auto const clockData = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(event);
625  compareSpacePointsToSim(event, clockData, recobHitToVoxelIDMap, voxelIDToPlaneTDCIDEMap);
626 
627  // Make sure the output tuples are filled
628  fHitSpacePointObj.fill();
629 
630  for(auto& hitObj : fHitSimObjVec) hitObj.fill();
631 
632  return;
633 }
634 
635 void SpacePointAnalysisMC::makeTrackToChanChargeMap(const TrackIDChanToTDCIDEMap& trackIDChanToTDCIDEMap,
636  TrackToChanChargeMap& trackToChanChargeMap,
637  float& bestTotDepEne,
638  int& bestTrackID) const
639 {
640  // Pretty straightforward looping here...
641  for(const auto& trackIDEPair : trackIDChanToTDCIDEMap)
642  {
643  ChanToChargeMap& chanToChargeMap = trackToChanChargeMap[trackIDEPair.first];
644 
645  float trackTotDepE(0.);
646 
647  for(const auto& chanTDCIDEPair : trackIDEPair.second)
648  {
649  ChargeDepositVec& chargeDepositVec = chanToChargeMap[chanTDCIDEPair.first];
650 
651  // Keep track of first,peak,last/ene
652  TDCIDEPair firstPair = chanTDCIDEPair.second.front();
653  TDCIDEPair peakPair = firstPair;
654  TDCIDEPair lastPair = chanTDCIDEPair.second.back();
655 
656  // Keep watch for gaps
657  TDCIDEPair prevPair = firstPair;
658 
659  // Keep track of deposited energy on a snippet
660  float snippetDepEne(0.);
661  float snippetNumElectrons(0.);
662 
663  for(const auto& tdcIDEPair : chanTDCIDEPair.second)
664  {
665  float depEne = tdcIDEPair.second->energy;
666 
667  trackTotDepE += depEne;
668 
669  // Watch for a gap...
670  if (tdcIDEPair.first - prevPair.first > 1)
671  {
672  chargeDepositVec.emplace_back(firstPair,peakPair,prevPair,snippetDepEne,snippetNumElectrons);
673 
674  firstPair = tdcIDEPair;
675  peakPair = firstPair;
676  snippetDepEne = 0.;
677  snippetNumElectrons = 0.;
678  }
679 
680  snippetDepEne += depEne;
681  snippetNumElectrons += tdcIDEPair.second->numElectrons;
682 
683  if (depEne > peakPair.second->energy) peakPair = tdcIDEPair;
684 
685  prevPair = tdcIDEPair;
686  }
687 
688  chargeDepositVec.emplace_back(firstPair,peakPair,lastPair,snippetDepEne,snippetNumElectrons);
689  }
690 
691  if (trackTotDepE > bestTotDepEne)
692  {
693  bestTrackID = trackIDEPair.first;
694  bestTotDepEne = trackTotDepE;
695  }
696  }
697 
698  return;
699 }
700 
701 void SpacePointAnalysisMC::compareHitsToSim(const art::Event& event, // For recovering data from event store
702  const ChanToTDCToIDEMap& chanToTDCToIDEMap, // This gives us ability to retrieve total charge deposits
703  const ChanToChargeMap& chanToChargeMap, // Charge deposit for specific track
704  const ChanToTDCIDEMap& chanToTDCIDEMap, // Charge deposit for specific track
705  const IDEToVoxelIDMap& ideToVoxelIDMap, // Mapping of ide info to voxels
706  RecobHitToVoxelIDMap& recobHitToVoxelIDMap) const // The output info
707 {
708  // We start by building a mapping between channels and lists of hits on that channel (ordered by time)
709  using HitPointerVec = std::vector<const recob::Hit*>;
710  using ChanToHitVecMap = std::map<raw::ChannelID_t,HitPointerVec>;
711 
712  ChanToHitVecMap chanToHitVecMap;
713 
714  // And now fill it
715  for(const auto& hitLabel : fRecobHitLabelVec)
716  {
717  art::Handle< std::vector<recob::Hit> > hitHandle;
718  event.getByLabel(hitLabel, hitHandle);
719 
720  // If no hits then skip
721  if ((*hitHandle).empty()) continue;
722 
723  for(const auto& hit : *hitHandle) chanToHitVecMap[hit.Channel()].emplace_back(&hit);
724  }
725 
726  // Now go through and order each vector of hits by time
727  for(auto& chanToHitPair : chanToHitVecMap)
728  {
729  HitPointerVec& hitPtrVec = chanToHitPair.second;
730 
731  std::sort(hitPtrVec.begin(),
732  hitPtrVec.end(),
733  [](const auto& left, const auto& right){return left->Channel() == right->Channel() ? left->PeakTime() < right->PeakTime() : left->Channel() < right->Channel();});
734  }
735 
736  auto const clockData = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(event);
737 
738  // The idea is to loop over the input sim information so we can look at efficiency as well as resolution issues
739  for(const auto& chanToChargePair : chanToChargeMap)
740  {
741  // Recover the channel
742  raw::ChannelID_t channel = chanToChargePair.first;
743 
744  // Look up the hits associated to this channel
745  ChanToHitVecMap::const_iterator chanToHitVecItr = chanToHitVecMap.find(channel);
746 
747  // For now we simply punt...
748  if (chanToHitVecItr == chanToHitVecMap.end()) continue;
749 
750  // Recover channel information based on this hit
751  const ChargeDepositVec& chargeDepositVec = chanToChargePair.second;
752 
753  // Get the hits...
754  const HitPointerVec& hitPtrVec = chanToHitVecItr->second;
755 
756  short int lastSnippetStart(-1);
757 
758  HitPointerVec hitVec;
759 
760  // Outer loop over hits in this hit collection
761  for(const auto& hitPtr : hitPtrVec)
762  {
763  // We want to collect together the hits that are on the same snippet. Hits will come grouped and in order
764  // along the snippet, so we simply keep them in a local vector until we hit the end of the snippet...
765  //
766  // ** It is worth noting this scheme as implemented will miss the last snippet of hits... so think about
767  // that for the future
768  short int snippetStart = hitPtr->StartTick();
769 
770  if (snippetStart == lastSnippetStart)
771  {
772  lastSnippetStart = snippetStart;
773 
774  hitVec.emplace_back(hitPtr);
775 
776  continue;
777  }
778 
779  // Process the current list of hits (which will be on the same snippet)
780  matchHitSim(clockData, hitVec, chanToTDCToIDEMap, chargeDepositVec, chanToTDCIDEMap, ideToVoxelIDMap, recobHitToVoxelIDMap);
781 
782  hitVec.clear();
783  hitVec.emplace_back(hitPtr);
784  lastSnippetStart = snippetStart;
785  }
786 
787  // Make sure to catch the last set of hits in the group
788  if (!hitVec.empty()) matchHitSim(clockData, hitVec, chanToTDCToIDEMap, chargeDepositVec, chanToTDCIDEMap, ideToVoxelIDMap, recobHitToVoxelIDMap);
789  }
790 
791  return;
792 }
793 
794 void SpacePointAnalysisMC::matchHitSim(const detinfo::DetectorClocksData& clockData,
795  const HitPointerVec& hitPointerVec, // Hits to match to simulation
796  const ChanToTDCToIDEMap& chanToTDCToIDEMap, // This gives us ability to retrieve total charge deposits
797  const ChargeDepositVec& chargeDepositVec, // Charge deposit for specific track
798  const ChanToTDCIDEMap& chanToTDCIDEMap, // Charge deposit for specific track
799  const IDEToVoxelIDMap& ideToVoxelIDMap, // Mapping of ide info to voxels
800  RecobHitToVoxelIDMap& recobHitToVoxelIDMap) const // The output info
801 {
802  // Data structure to allow ordering of multiple hits in a snippet
803  using HitPeakTimeChargeTuple = std::tuple<int,const recob::Hit*,ChargeDepositVec::const_iterator>;
804  using HitPeakTimeChargeVec = std::vector<HitPeakTimeChargeTuple>;
805 
806  HitPeakTimeChargeVec hitPeakTimeChargeVec;
807  float maxPulseHeight(1.); // Don't let this be zero
808 
809  // If here then we are on to the next hit, so we need to process our current list
810  for(const auto& hit : hitPointerVec)
811  {
812  // Recover hit time range (in ticks), cast a wide net here
813  int peakTick = std::round(hit->PeakTime());
814  int startTick = std::max( 0,int(std::floor(hit->PeakTime() - 3. * hit->RMS())));
815  int endTick = std::min(4096,int(std::ceil(hit->PeakTime() + 3. * hit->RMS())));
816 
817  int startTDC = clockData.TPCTick2TDC(startTick - fOffsetVec[hit->WireID().Plane]);
818  int peakTDC = clockData.TPCTick2TDC(peakTick - fOffsetVec[hit->WireID().Plane]);
819  int endTDC = clockData.TPCTick2TDC(endTick - fOffsetVec[hit->WireID().Plane]);
820 
821  maxPulseHeight = std::max(maxPulseHeight,hit->PeakAmplitude());
822 
823  // If we have a match then this iterator gets set to the matching values
824  ChargeDepositVec::const_iterator chargeMatchItr = chargeDepositVec.end();
825 
826  int bestPeakDiff = std::numeric_limits<int>::max();
827 
828  // Match the hit (if there is one)
829  for(ChargeDepositVec::const_iterator chargeInfoItr = chargeDepositVec.begin(); chargeInfoItr != chargeDepositVec.end(); chargeInfoItr++)
830  {
831  // Require some amount of overlap between the hit and the sim info
832  if (endTDC > std::get<0>(*chargeInfoItr).first && startTDC < std::get<2>(*chargeInfoItr).first)
833  {
834  const TDCIDEPair& peakTDCIDE = std::get<1>(*chargeInfoItr);
835 
836  int peakDiff = peakTDC - int(peakTDCIDE.first);
837 
838  if (std::abs(peakDiff) < std::abs(bestPeakDiff))
839  {
840  bestPeakDiff = peakDiff;
841  chargeMatchItr = chargeInfoItr;
842  }
843  }
844  }
845 
846  // If no match then skip
847  if (chargeMatchItr == chargeDepositVec.end()) continue;
848 
849  hitPeakTimeChargeVec.emplace_back(std::make_tuple(bestPeakDiff,hit,chargeMatchItr));
850  }
851 
852  if (!hitPeakTimeChargeVec.empty())
853  {
854  // Ok, now we sort this vector by smallest peak time
855  std::sort(hitPeakTimeChargeVec.begin(),hitPeakTimeChargeVec.end(),[](const auto& left,const auto& right){return std::abs(std::get<0>(left)) < std::abs(std::get<0>(right));});
856 
857  // Keep track of hit ordering on this snippet
858  int hitOrder(0);
859 
860  HitSimulationTupleObj& hitObj = fHitSimObjVec[std::get<1>(hitPeakTimeChargeVec.front())->WireID().Plane];
861 
862  // Now loop through
863  for(const auto& hitPeakCharge : hitPeakTimeChargeVec)
864  {
865  const recob::Hit* hit = std::get<1>(hitPeakCharge);
866  const ChargeDeposit& chargeDeposit = *std::get<2>(hitPeakCharge);
867 
868  // Recover hit time range (in ticks), cast a wide net here
869  int peakTick = std::round(hit->PeakTime());
870  int startTick = std::max( 0,int(std::floor(hit->PeakTime() - 3. * hit->RMS())));
871  int endTick = std::min(4096,int(std::ceil(hit->PeakTime() + 3. * hit->RMS())));
872 
873  int startTDC = clockData.TPCTick2TDC(startTick - fOffsetVec[hit->WireID().Plane]);
874  int peakTDC = clockData.TPCTick2TDC(peakTick - fOffsetVec[hit->WireID().Plane]);
875  int endTDC = clockData.TPCTick2TDC(endTick - fOffsetVec[hit->WireID().Plane]);
876 
877  int firstSimTick(std::get<0>(chargeDeposit).first);
878  int lastSimTick(std::get<2>(chargeDeposit).first);
879  int maxDepTick(std::get<1>(chargeDeposit).first);
880  float maxDepEneTick(std::get<1>(chargeDeposit).second->energy);
881  float bestNumElectrons(std::get<4>(chargeDeposit));
882  float bestDepEne(std::get<3>(chargeDeposit));
883  float totDepEne(0.);
884  float totNumElectrons(0.);
885  int bestTicks(lastSimTick - firstSimTick + 1);
886 
887  // We want to get the total energy deposit from all particles in the ticks for this hit
888  const TDCToIDEMap& tdcToIDEMap = chanToTDCToIDEMap.find(hit->Channel())->second;
889  for(const auto& tdcToIDEPair : tdcToIDEMap)
890  {
891  for(const auto& ide : tdcToIDEPair.second)
892  {
893  totDepEne += ide->energy;
894  totNumElectrons += ide->numElectrons;
895  }
896  }
897 
898  // One final time through to find sim ticks that "matter"
899  // We define this as the collection of IDE's that make up to 90% of the total deposit
900  ChanToTDCIDEMap::const_iterator tickToTDCIDEVecItr = chanToTDCIDEMap.find(hit->Channel());
901 
902  if (tickToTDCIDEVecItr == chanToTDCIDEMap.end()) continue;
903 
904  const TickTDCIDEVec& tickToTDCIDEVec = tickToTDCIDEVecItr->second;
905  TickTDCIDEVec tickIDEVec;
906 
907  for(const auto& tickInfo : tickToTDCIDEVec)
908  {
909  //if (tickInfo.first >= firstSimTick && tickInfo.first <= lastSimTick) tickIDEVec.emplace_back(tickInfo);
910  if (tickInfo.first >= startTDC && tickInfo.first <= endTDC) tickIDEVec.emplace_back(tickInfo);
911  }
912 
913  std::sort(tickIDEVec.begin(),tickIDEVec.end(),[](const auto& left,const auto& right){return left.second->energy > right.second->energy;});
914 
915  // Grab the voxelID set for this tick
916  VoxelIDSet voxelIDSet;
917 
918  int bestTicksGood(0);
919  float sumEne(0.);
920 
921  for(const auto& tickInfo : tickIDEVec)
922  {
923  // At the same time we can keep track of the voxels associated to the best track
924  IDEToVoxelIDMap::const_iterator ideToVoxelIDMapItr = ideToVoxelIDMap.find(tickInfo.second);
925 
926  if (ideToVoxelIDMapItr == ideToVoxelIDMap.end()) continue;
927 
928  const sim::LArVoxelID& voxelID = ideToVoxelIDMapItr->second;
929 
930  sumEne += tickInfo.second->energy;
931  bestTicksGood++;
932 
933  voxelIDSet.insert(voxelID);
934 
935  if (sumEne > 0.9 * bestDepEne) break;
936  }
937 
938  // Finally, grab the voxels from the track leaving the most energy
939  recobHitToVoxelIDMap[hit] = voxelIDSet;
940 
941  hitObj.fillSimInfo(bestTicks, bestTicksGood, totDepEne, totNumElectrons, bestDepEne, bestNumElectrons, maxDepEneTick);
942  hitObj.fillMixedInfo(endTick-startTick+1, maxDepTick-startTDC, peakTDC-maxDepTick);
943  hitObj.fillHitInfo(hit,hitOrder++);
944  }
945 
946  // Resort in pulse height order (largest to smallest)
947  std::sort(hitPeakTimeChargeVec.begin(),hitPeakTimeChargeVec.end(),[](const auto& left,const auto& right){return std::get<1>(left)->PeakAmplitude() > std::get<1>(right)->PeakAmplitude();});
948 
949  // Now loop through
950  for(const auto& hitPeakCharge : hitPeakTimeChargeVec)
951  {
952  hitObj.fPHOrderHitVec.emplace_back(std::get<1>(hitPeakCharge)->LocalIndex());
953  hitObj.fPHFractionVec.emplace_back(std::get<1>(hitPeakCharge)->PeakAmplitude() / maxPulseHeight);
954  }
955  }
956 
957  return;
958 }
959 
960 void SpacePointAnalysisMC::compareSpacePointsToSim(const art::Event& event,
961  const detinfo::DetectorClocksData& clockData,
962  const RecobHitToVoxelIDMap& recobHitToVoxelIDMap,
963  const VoxelIDToPlaneTDCIDEMap& voxelToPlaneTDCIDEMap) const
964 {
965  // Armed with these maps we can now process the SpacePoints...
966  if (!recobHitToVoxelIDMap.empty())
967  {
968  // Diagnostics
969  using Triplet = std::tuple<const recob::Hit*, const recob::Hit*, const recob::Hit*>;
970  using TripletMap = std::map<Triplet,std::vector<const recob::SpacePoint*>>;
971 
972  TripletMap tripletMap;
973 
974  // So now we loop through the various SpacePoint sources
975  for(const auto& spacePointLabel : fSpacePointLabelVec)
976  {
977  art::Handle< std::vector<recob::SpacePoint>> spacePointHandle;
978  event.getByLabel(spacePointLabel, spacePointHandle);
979 
980  if (!spacePointHandle.isValid()) continue;
981 
982  // Look up assocations to hits
983  art::FindManyP<recob::Hit> spHitAssnVec(spacePointHandle, event, spacePointLabel);
984 
985  // And now, without further adieu, here we begin the loop that will actually produce some useful output.
986  // Loop over all space points and find out their true nature
987  for(size_t idx = 0; idx < spacePointHandle->size(); idx++)
988  {
989  art::Ptr<recob::SpacePoint> spacePointPtr(spacePointHandle,idx);
990 
991  std::vector<art::Ptr<recob::Hit>> associatedHits(spHitAssnVec.at(spacePointPtr.key()));
992 
993  if (associatedHits.size() != 3)
994  {
995  mf::LogDebug("SpacePointAnalysisMC") << "I am certain this cannot happen... but here you go, space point with " << associatedHits.size() << " hits" << std::endl;
996  continue;
997  }
998 
999  // Retrieve the magic numbers from the space point
1000  float spQuality = spacePointPtr->Chisq();
1001  float spCharge = spacePointPtr->ErrXYZ()[1];
1002  float spAsymmetry = spacePointPtr->ErrXYZ()[3];
1003  float smallestPH = std::numeric_limits<float>::max();
1004  float largestPH = 0.;
1005  int numHits = 0;
1006  float averagePH = 0.;
1007  float averagePT = 0.;
1008  float largestDelT = 0.;
1009 
1010  std::vector<int> numIDEsHitVec;
1011  std::vector<int> simHitDeltaTVec = {0,0,0};
1012  std::vector<float> hitPeakTimeVec = {-100.,-100.,-100.};
1013  std::vector<float> bigElecDepVec = {0.,0.,0.};
1014  std::vector<unsigned short> bigTDCVec = {0,0,0};
1015  int numIDEsSpacePoint(0);
1016  int numLongHits(0);
1017  int numIntersections(0);
1018 
1019  std::vector<RecobHitToVoxelIDMap::const_iterator> recobHitToVoxelIterVec;
1020 
1021  std::vector<const recob::Hit*> recobHitVec(3,nullptr);
1022 
1023  // Now we can use our maps to find out if the hits making up the SpacePoint are truly related...
1024  for(const auto& hitPtr : associatedHits)
1025  {
1026  RecobHitToVoxelIDMap::const_iterator hitToVoxelItr = recobHitToVoxelIDMap.find(hitPtr.get());
1027 
1028  float peakAmplitude = hitPtr->PeakAmplitude();
1029  float peakTime = hitPtr->PeakTime();
1030  int plane = hitPtr->WireID().Plane;
1031 
1032  recobHitVec[plane] = hitPtr.get();
1033 
1034  numHits++;
1035  averagePH += peakAmplitude;
1036  averagePT += peakTime;
1037 
1038  smallestPH = std::min(peakAmplitude,smallestPH);
1039  largestPH = std::max(peakAmplitude,largestPH);
1040 
1041  if (hitPtr->DegreesOfFreedom() < 2) numLongHits++;
1042 
1043  if (hitToVoxelItr == recobHitToVoxelIDMap.end())
1044  {
1045  numIDEsHitVec.push_back(0);
1046  continue;
1047  }
1048 
1049  recobHitToVoxelIterVec.push_back(hitToVoxelItr);
1050  numIDEsHitVec.push_back(hitToVoxelItr->second.size());
1051 
1052  hitPeakTimeVec[plane] = clockData.TPCTick2TDC(peakTime);
1053  }
1054 
1055  Triplet hitTriplet(recobHitVec[0],recobHitVec[1],recobHitVec[2]);
1056 
1057  tripletMap[hitTriplet].emplace_back(spacePointPtr.get());
1058 
1059  averagePH /= float(numHits);
1060  averagePT /= float(numHits);
1061 
1062  for(const auto& hitPtr : associatedHits)
1063  {
1064  float delT = hitPtr->PeakTime() - averagePT;
1065 
1066  if (std::abs(delT) > std::abs(largestDelT)) largestDelT = delT;
1067  }
1068 
1069  // If a SpacePoint is made from "true" MC hits then we will have found the relations to the MC info for all three
1070  // hits. If this condition is not satisfied it means one or more hits making the SpacePoint are noise hits
1071  if (recobHitToVoxelIterVec.size() == 3)
1072  {
1073  // Find the intersection of the vectors of IDEs for the first two hits
1074  std::vector<sim::LArVoxelID> firstIntersectionVec(recobHitToVoxelIterVec[0]->second.size()+recobHitToVoxelIterVec[1]->second.size());
1075 
1076  std::vector<sim::LArVoxelID>::iterator firstIntersectionItr = std::set_intersection(recobHitToVoxelIterVec[0]->second.begin(),recobHitToVoxelIterVec[0]->second.end(),
1077  recobHitToVoxelIterVec[1]->second.begin(),recobHitToVoxelIterVec[1]->second.end(),
1078  firstIntersectionVec.begin());
1079 
1080  firstIntersectionVec.resize(firstIntersectionItr - firstIntersectionVec.begin());
1081 
1082  // No intersection means, of course, the hits did not come from the same MC energy deposit
1083  if (!firstIntersectionVec.empty())
1084  {
1085  // Now find the intersection of the resultant intersection above and the third hit
1086  std::vector<sim::LArVoxelID> secondIntersectionVec(firstIntersectionVec.size()+recobHitToVoxelIterVec[2]->second.size());
1087 
1088  std::vector<sim::LArVoxelID>::iterator secondIntersectionItr = std::set_intersection(firstIntersectionVec.begin(), firstIntersectionVec.end(),
1089  recobHitToVoxelIterVec[2]->second.begin(),recobHitToVoxelIterVec[2]->second.end(),
1090  secondIntersectionVec.begin());
1091 
1092  secondIntersectionVec.resize(secondIntersectionItr - secondIntersectionVec.begin());
1093 
1094  numIntersections++;
1095 
1096  // Again, no IDEs in the intersection means it is a ghost space point but, of course, we are hoping
1097  // there are common IDEs so we can call it a real SpacePoint
1098  if (!secondIntersectionVec.empty())
1099  {
1100  for(const sim::LArVoxelID& voxelID : secondIntersectionVec)
1101  {
1102  VoxelIDToPlaneTDCIDEMap::const_iterator planeToTDCToIDESetMap = voxelToPlaneTDCIDEMap.find(voxelID);
1103 
1104  if (planeToTDCToIDESetMap->second.size() > 2)
1105  {
1106  numIDEsSpacePoint += 1;
1107 
1108  for(const auto& planeInfoPair : planeToTDCToIDESetMap->second)
1109  {
1110  unsigned short plane = planeInfoPair.first;
1111  float phBig = 0.;
1112  unsigned short tdcBig = 0;
1113 
1114  for(const auto& tdcIDEPair : planeInfoPair.second)
1115  {
1116  for(const auto& ide : tdcIDEPair.second)
1117  {
1118  if (phBig < ide->numElectrons)
1119  {
1120  phBig = ide->numElectrons;
1121  tdcBig = tdcIDEPair.first;
1122  }
1123  }
1124  }
1125 
1126  bigElecDepVec[plane] = phBig;
1127  bigTDCVec[plane] = tdcBig;
1128  }
1129  }
1130  else std::cout << " --> Not matching all three planes" << std::endl;
1131  }
1132 
1133  numIntersections++;
1134  }
1135  }
1136  }
1137 
1138  // Fill for "all" cases
1139  fHitSpacePointObj.fSPQualityVec.push_back(spQuality);
1140  fHitSpacePointObj.fSPTotalChargeVec.push_back(spCharge);
1141  fHitSpacePointObj.fSPAsymmetryVec.push_back(spAsymmetry);
1142  fHitSpacePointObj.fSmallestPHVec.push_back(smallestPH);
1143  fHitSpacePointObj.fLargestPHVec.push_back(largestPH);
1144  fHitSpacePointObj.fAveragePHVec.push_back(averagePH);
1145  fHitSpacePointObj.fLargestDelTVec.push_back(largestDelT);
1146 
1147  fHitSpacePointObj.fNumIDEsHit0Vec.push_back(numIDEsHitVec[0]);
1148  fHitSpacePointObj.fNumIDEsHit1Vec.push_back(numIDEsHitVec[1]);
1149  fHitSpacePointObj.fNumIDEsHit2Vec.push_back(numIDEsHitVec[2]);
1150  fHitSpacePointObj.fNumIDEsSpacePointVec.push_back(numIDEsSpacePoint);
1151 
1152  fHitSpacePointObj.fNumLongHitsVec.emplace_back(numLongHits);
1153  fHitSpacePointObj.fNumPlanesSimMatchVec.emplace_back(recobHitToVoxelIterVec.size());
1154  fHitSpacePointObj.fNumIntersectSetVec.emplace_back(numIntersections);
1155  fHitSpacePointObj.fSimHitDeltaT0Vec.emplace_back(bigTDCVec[0] - hitPeakTimeVec[0]);
1156  fHitSpacePointObj.fSimHitDeltaT1Vec.emplace_back(bigTDCVec[1] - hitPeakTimeVec[1]);
1157  fHitSpacePointObj.fSimHitDeltaT2Vec.emplace_back(bigTDCVec[2] - hitPeakTimeVec[2]);
1158 
1159  fHitSpacePointObj.fSimDelta10Vec.emplace_back(bigTDCVec[1] - bigTDCVec[0]);
1160  fHitSpacePointObj.fSimDelta21Vec.emplace_back(bigTDCVec[2] - bigTDCVec[1]);
1161  fHitSpacePointObj.fHitDelta10Vec.emplace_back(hitPeakTimeVec[1] - hitPeakTimeVec[0]);
1162  fHitSpacePointObj.fHitDelta21Vec.emplace_back(hitPeakTimeVec[2] - hitPeakTimeVec[1]);
1163  fHitSpacePointObj.fBigElecDep0Vec.emplace_back(bigElecDepVec[0]);
1164  fHitSpacePointObj.fBigElecDep1Vec.emplace_back(bigElecDepVec[1]);
1165  fHitSpacePointObj.fBigElecDep2Vec.emplace_back(bigElecDepVec[2]);
1166  }
1167  }
1168 
1169  // Can we check to see if we have duplicates?
1170  std::vector<int> numSpacePointVec = {0,0,0,0,0};
1171  for(const auto& tripletPair : tripletMap)
1172  {
1173  int numSpacePoints = std::min(numSpacePointVec.size()-1,tripletPair.second.size());
1174  numSpacePointVec[numSpacePoints]++;
1175  }
1176  std::cout << "====>> Found " << tripletMap.size() << " SpacePoints, numbers: ";
1177  for(const auto& count : numSpacePointVec) std::cout << count << " ";
1178  std::cout << std::endl;
1179  }
1180 
1181  return;
1182 }
1183 
1184 
1185 // Useful for normalizing histograms
1186 void SpacePointAnalysisMC::endJob(int numEvents)
1187 {
1188  return;
1189 }
1190 
1191 DEFINE_ART_CLASS_TOOL(SpacePointAnalysisMC)
1192 }
recob::Hit::LocalIndex
short int LocalIndex() const
How well do we believe we know this hit?
Definition: Hit.h:227
SpacePointAnalysisMC::SpacePointAnalysisMC::endJob
void endJob(int numEvents) override
Interface for method to executve at the end of run processing.
Definition: SpacePointAnalysisMC_tool.cc:1186
SpacePointAnalysisMC::SpacePointAnalysisMC::TrackToChanChargeMap
std::unordered_map< int, ChanToChargeMap > TrackToChanChargeMap
Definition: SpacePointAnalysisMC_tool.cc:390
SpacePointAnalysisMC::HitSimulationTupleObj::fillMixedInfo
void fillMixedInfo(int hitWidth, int ticksToMax, int deltaTicks)
Definition: SpacePointAnalysisMC_tool.cc:260
second
second
Definition: detectorproperties_icarus.fcl:8
ServiceUtil.h
Utilities related to art service access.
SpacePointAnalysisMC::SpacePointAnalysisMC::PlaneToTDCToIDESetMap
std::map< unsigned short, TDCToIDESetMap > PlaneToTDCToIDESetMap
Definition: SpacePointAnalysisMC_tool.cc:377
SpacePointAnalysisMC::SpacePointAnalysisMC::TrackIDChanToTDCIDEMap
std::unordered_map< int, ChanToTDCIDEMap > TrackIDChanToTDCIDEMap
Definition: SpacePointAnalysisMC_tool.cc:384
recob::Hit::WireID
geo::WireID WireID() const
Definition: Hit.h:233
recob::Hit::RMS
float RMS() const
RMS of the hit shape, in tick units.
Definition: Hit.h:220
SpacePointAnalysisMC::SpacePointAnalysisMC::fillHistograms
void fillHistograms(const art::Event &) const override
Interface for filling histograms.
Definition: SpacePointAnalysisMC_tool.cc:533
SpacePointAnalysisMC::SpacePointAnalysisMC::initializeTuple
void initializeTuple(TTree *) override
Interface for initializing the tuple variables.
Definition: SpacePointAnalysisMC_tool.cc:488
SpacePointAnalysisMC::HitSimulationTupleObj::fillSimInfo
void fillSimInfo(int ticksSimChannel, int ticksSimChanMost, float totDepEne, float totNumElectrons, float bestDepEne, float bestNumElectrons, float maxDepEneTick)
Definition: SpacePointAnalysisMC_tool.cc:242
SpacePointAnalysisMC::SpacePointAnalysisMC::ChargeDeposit
std::tuple< TDCIDEPair, TDCIDEPair, TDCIDEPair, float, float > ChargeDeposit
Definition: SpacePointAnalysisMC_tool.cc:387
Hit.h
Declaration of signal hit object.
SpacePointAnalysisMC::SpacePointAnalysisMC::fOffsetVec
std::vector< int > fOffsetVec
Allow offsets for each plane.
Definition: SpacePointAnalysisMC_tool.cc:415
right
walls no right
Definition: selectors.fcl:105
SpacePointAnalysisMC::HitSimulationTupleObj::fillHitInfo
void fillHitInfo(const recob::Hit *hit, int hitOrder)
Definition: SpacePointAnalysisMC_tool.cc:269
SpacePointAnalysisMC::SpacePointAnalysisMC::matchHitSim
void matchHitSim(const detinfo::DetectorClocksData &clockData, const HitPointerVec &, const ChanToTDCToIDEMap &, const ChargeDepositVec &, const ChanToTDCIDEMap &, const IDEToVoxelIDMap &, RecobHitToVoxelIDMap &) const
Definition: SpacePointAnalysisMC_tool.cc:794
recob::Hit::DegreesOfFreedom
int DegreesOfFreedom() const
Definition: Hit.h:229
geo::GeometryCore::ChannelToWire
std::vector< geo::WireID > ChannelToWire(raw::ChannelID_t const channel) const
Returns a list of wires connected to the specified TPC channel.
Definition: GeometryCore.cxx:1196
SpacePointAnalysisMC::SpacePointAnalysisMC::TDCIDEPair
std::pair< unsigned short, const sim::IDE * > TDCIDEPair
Definition: SpacePointAnalysisMC_tool.cc:381
RawDigit.h
Definition of basic raw digits.
recob::Hit::Integral
float Integral() const
Integral under the calibrated signal waveform of the hit, in tick x ADC units.
Definition: Hit.h:224
hit
process_name hit
Definition: cheaterreco.fcl:51
recob::Hit::GoodnessOfFit
float GoodnessOfFit() const
Degrees of freedom in the determination of the hit signal shape (-1 by default)
Definition: Hit.h:228
recob::Hit::Multiplicity
short int Multiplicity() const
How many hits could this one be shared with.
Definition: Hit.h:226
SpacePointAnalysisMC::SpacePointAnalysisMC::compareHitsToSim
void compareHitsToSim(const art::Event &, const ChanToTDCToIDEMap &, const ChanToChargeMap &, const ChanToTDCIDEMap &, const IDEToVoxelIDMap &, RecobHitToVoxelIDMap &) const
Definition: SpacePointAnalysisMC_tool.cc:701
recob::Hit::PeakAmplitude
float PeakAmplitude() const
The estimated amplitude of the hit at its peak, in ADC units.
Definition: Hit.h:221
SpacePointAnalysisMC::SpacePointAnalysisMC::initializeHists
void initializeHists(art::ServiceHandle< art::TFileService > &, const std::string &) override
Interface for initializing the histograms to be filled.
Definition: SpacePointAnalysisMC_tool.cc:480
abs
T abs(T value)
Definition: sparse_vector_test.cc:30
SpacePointAnalysisMC::SpacePointAnalysisMC::TDCToIDEMap
std::map< unsigned short, SimIDESet > TDCToIDEMap
Definition: SpacePointAnalysisMC_tool.cc:371
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
geo::fhicl::WireID
IDparameter< geo::WireID > WireID
Member type of validated geo::WireID parameter.
Definition: geo_types_fhicl.h:407
SpacePointAnalysisMC::SpacePointAnalysisMC::SimIDESet
std::set< const sim::IDE *, ideCompare > SimIDESet
Definition: SpacePointAnalysisMC_tool.cc:368
sim::IDE
Ionization at a point of the TPC sensitive volume.
Definition: SimChannel.h:86
SpacePointAnalysisMC::SpacePointAnalysisMC::VoxelIDToPlaneTDCIDEMap
std::map< sim::LArVoxelID, PlaneToTDCToIDESetMap > VoxelIDToPlaneTDCIDEMap
Definition: SpacePointAnalysisMC_tool.cc:378
SpacePointAnalysisMC::SpacePointAnalysisMC::compareSpacePointsToSim
void compareSpacePointsToSim(const art::Event &, const detinfo::DetectorClocksData &clockData, const RecobHitToVoxelIDMap &, const VoxelIDToPlaneTDCIDEMap &) const
Definition: SpacePointAnalysisMC_tool.cc:960
sim::IDE::energy
float energy
energy deposited by ionization by this track ID and time [MeV]
Definition: SimChannel.h:118
SpacePointAnalysisMC::SpacePointAnalysisMC::ChanToChargeMap
std::map< raw::ChannelID_t, ChargeDepositVec > ChanToChargeMap
Definition: SpacePointAnalysisMC_tool.cc:389
SpacePointAnalysisMC::SpacePointAnalysisMC::fGeometry
const geo::GeometryCore * fGeometry
pointer to Geometry service
Definition: SpacePointAnalysisMC_tool.cc:432
SpacePointAnalysisMC::SpacePointAnalysisMC::HitSimObjVec
std::vector< HitSimulationTupleObj > HitSimObjVec
Definition: SpacePointAnalysisMC_tool.cc:426
SpacePointAnalysisMC::SpacePointAnalysisMC::ChanToTDCToIDEMap
std::map< raw::ChannelID_t, TDCToIDEMap > ChanToTDCToIDEMap
Definition: SpacePointAnalysisMC_tool.cc:372
left
walls no left
Definition: selectors.fcl:105
SpacePointAnalysisMC::HitPtrVec
std::vector< art::Ptr< recob::Hit >> HitPtrVec
Definition: SpacePointAnalysisMC_tool.cc:62
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
SpacePointAnalysisMC::SpacePointAnalysisMC::ideCompare::operator()
bool operator()(const sim::IDE *left, const sim::IDE *right) const
Definition: SpacePointAnalysisMC_tool.cc:364
SpacePointAnalysisMC::SpacePointAnalysisMC::IDEToVoxelIDMap
std::unordered_map< const sim::IDE *, sim::LArVoxelID > IDEToVoxelIDMap
Definition: SpacePointAnalysisMC_tool.cc:369
recob::Hit::PeakTime
float PeakTime() const
Time of the signal peak, in tick units.
Definition: Hit.h:218
SpacePointAnalysisMC::TrackViewHitMap
std::map< int, ViewHitMap > TrackViewHitMap
Definition: SpacePointAnalysisMC_tool.cc:64
SpacePointAnalysisMC::ViewHitMap
std::map< size_t, HitPtrVec > ViewHitMap
Definition: SpacePointAnalysisMC_tool.cc:63
detinfo::DetectorClocksData
Contains all timing reference information for the detector.
Definition: DetectorClocksData.h:283
icarus::trigger::to_string
std::string to_string(WindowPattern const &pattern)
Definition: WindowPattern.h:115
SpacePointAnalysisMC::SpacePointAnalysisMC::configure
void configure(fhicl::ParameterSet const &pset) override
Definition: SpacePointAnalysisMC_tool.cc:464
SimEnergyDeposit.h
contains information for a single step in the detector simulation
SpacePointAnalysisMC::SpacePointAnalysisMC::TDCToIDESetMap
std::unordered_map< unsigned short, SimIDESet > TDCToIDESetMap
Definition: SpacePointAnalysisMC_tool.cc:376
SimChannel.h
object containing MC truth information necessary for making RawDigits and doing back tracking ...
ChannelStatusProvider.h
Interface for experiment-specific channel quality info provider.
SpacePointAnalysisMC::SpacePointAnalysisMC::VoxelIDToIDESetMap
std::map< sim::LArVoxelID, SimIDESet > VoxelIDToIDESetMap
Definition: SpacePointAnalysisMC_tool.cc:370
recob::Hit::SummedADC
float SummedADC() const
The sum of calibrated ADC counts of the hit (0. by default)
Definition: Hit.h:223
Wire.h
Declaration of basic channel signal object.
SpacePointAnalysisMC::SpacePointAnalysisMC::ChanToTDCIDEMap
std::unordered_map< raw::ChannelID_t, TickTDCIDEVec > ChanToTDCIDEMap
Definition: SpacePointAnalysisMC_tool.cc:383
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
count
std::size_t count(Cont const &cont)
Definition: IntegerRanges_test.cc:29
SpacePointAnalysisMC::SpacePointAnalysisMC::SpacePointAnalysisMC
SpacePointAnalysisMC(fhicl::ParameterSet const &pset)
Constructor.
Definition: SpacePointAnalysisMC_tool.cc:442
raw::ChannelID_t
unsigned int ChannelID_t
Type representing the ID of a readout channel.
Definition: RawTypes.h:28
ChannelStatusService.h
Interface for experiment-specific service for channel quality info.
LArVoxelID.h
Unique identifier for a given LAr voxel.
SpacePointAnalysisMC::SpacePointAnalysisMC::makeTrackToChanChargeMap
void makeTrackToChanChargeMap(const TrackIDChanToTDCIDEMap &, TrackToChanChargeMap &, float &, int &) const
Definition: SpacePointAnalysisMC_tool.cc:635
recob::Hit::Channel
raw::ChannelID_t Channel() const
ID of the readout channel the hit was extracted from.
Definition: Hit.h:230
detinfo::DetectorClocksData::TPCTick2TDC
double TPCTick2TDC(double const tick) const
Definition: DetectorClocksData.h:493
SpacePointAnalysisMC::SpacePointAnalysisMC::RecobHitToVoxelIDMap
std::unordered_map< const recob::Hit *, VoxelIDSet > RecobHitToVoxelIDMap
Definition: SpacePointAnalysisMC_tool.cc:397
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