LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
MCHitAnaExample_module.cc
Go to the documentation of this file.
1 // Class: MCHitAnaExample
3 // Module Type: analyzer
4 // File: MCHitAnaExample_module.cc
5 //
6 // Generated at Tue Jun 24 22:54:57 2014 by Kazuhiro Terao using artmod
7 // from cetpkgsupport v1_05_04.
9 
15 #include "art_root_io/TFileService.h"
16 #include "fhiclcpp/ParameterSet.h"
17 
22 #include <TH1D.h>
23 #include <TH2D.h>
24 #include <TStopwatch.h>
25 
26 namespace hit {
27 
28  class MCHitAnaExample;
29 
31  public:
32  explicit MCHitAnaExample(fhicl::ParameterSet const& p);
33  virtual ~MCHitAnaExample();
34 
35  private:
36  void analyze(art::Event const& e) override;
37 
38  //
39  // Module names
40  //
41 
43  std::string fRecoHitModuleName;
44 
46  std::string fMCHitModuleName;
47 
48  //
49  // Stopwatch for timings
50  //
52 
53  //
54  // Histogram for timing record
55  //
56 
59 
62 
65 
68 
71 
72  //
73  // Histograms purely from MCHit
74  //
75 
77  std::vector<TH1D*> hMCHitQ_v;
78 
80  std::vector<TH1D*> hMCHitMult_v;
81 
82  //
83  // Histograms purely from RecoHit
84  //
85 
87  std::vector<TH1D*> hRecoHitQ_v;
88 
90  std::vector<TH1D*> hRecoHitMult_v;
91 
92  //
93  // Histograms from MC+Reco combined/compared information
94  //
95 
97  std::vector<TH2D*> hCorrMult_v;
98 
100  std::vector<TH1D*> hVoidHitQ_v;
101 
103  std::vector<TH2D*> hCorrQ_v;
104 
106  std::vector<TH2D*> hCorrQSum_v;
107 
109  std::vector<TH1D*> hQRatio_v;
110 
115  std::vector<TH1D*> hQSumRatio_v;
116 
118  std::vector<TH1D*> hDT_v;
119 
121  std::vector<TH1D*> hMCHitMultPerRecoHit_v;
122  };
123 
125 
126  {
127  fRecoHitModuleName = p.get<std::string>("RecoHitModuleName");
128  fMCHitModuleName = p.get<std::string>("MCHitModuleName");
129 
131 
133 
134  for (unsigned char plane = 0; plane < geo->Nplanes(); ++plane) {
135 
136  hMCHitQ_v.push_back(
137  fs->make<TH1D>(Form("hMCHitQ_%d", plane),
138  Form("MCHit Charge (Plane %d); Charge [ADC]; # MCHit", plane),
139  250,
140  0,
141  2000));
142 
143  hMCHitMult_v.push_back(
144  fs->make<TH1D>(Form("hMCHitMult_%d", plane),
145  Form("MCHit Multiplicity (Plane %d); Multiplicity; # MCHit", plane),
146  250,
147  0,
148  5000));
149 
150  hRecoHitQ_v.push_back(
151  fs->make<TH1D>(Form("hRecoHitQ_%d", plane),
152  Form("RecoHit Charge (Plane %d); Charge [ADC]; # RecoHits", plane),
153  250,
154  0,
155  2000));
156 
157  hRecoHitMult_v.push_back(
158  fs->make<TH1D>(Form("hRecoHitMult_%d", plane),
159  Form("RecoHit Multiplicity (Plane %d); Multiplicity; # RecoHits", plane),
160  200,
161  0,
162  2000));
163 
164  hCorrMult_v.push_back(
165  fs->make<TH2D>(Form("hCorrMult_%d", plane),
166  Form("# MCHits vs. # RecoHits (Plane %d); # MCHits; # RecoHits", plane),
167  250,
168  0,
169  5000,
170  200,
171  0,
172  2000));
173 
174  hVoidHitQ_v.push_back(fs->make<TH1D>(
175  Form("hVoidHitQ_%d", plane),
176  Form("RecoHit Charge (No Corresponding MCHit)) (Plane %d); Charge [ADC]; # RecoHits",
177  plane),
178  250,
179  0,
180  2000));
181 
182  hCorrQ_v.push_back(fs->make<TH2D>(
183  Form("hCorrQ_%d", plane),
184  Form("RecoHit Q vs. Closest MCHit Q (Plane %d); MCHit Charge [ADC]; RecoHit Charge [ADC]",
185  plane),
186  250,
187  0,
188  2000,
189  250,
190  0,
191  2000));
192 
193  hCorrQSum_v.push_back(
194  fs->make<TH2D>(Form("hCorrQSum_%d", plane),
195  Form("RecoHit Q vs. MCHit QSum Within Start=>End (Plane %d); MCHit Charge "
196  "[ADC]; RecoHit Charge [ADC]",
197  plane),
198  250,
199  0,
200  2000,
201  250,
202  0,
203  2000));
204 
205  hQRatio_v.push_back(
206  fs->make<TH1D>(Form("hQRatio_%d", plane),
207  Form("Reco/MCHit Charge (Plane %d); Ratio; # RecoHits", plane),
208  200,
209  0,
210  2));
211 
212  hQSumRatio_v.push_back(
213  fs->make<TH1D>(Form("hQSumRatio_%d", plane),
214  Form("Reco/MCHit Charge (Plane %d); Ratio; # RecoHits", plane),
215  200,
216  0,
217  2));
218 
219  hDT_v.push_back(fs->make<TH1D>(
220  Form("hDT_%d", plane),
221  Form("#DeltaT btw Reco and Closest MCHit (Plane %d); #DeltaT; # RecoHits", plane),
222  190,
223  -9.5,
224  9.5));
225 
226  hMCHitMultPerRecoHit_v.push_back(fs->make<TH1D>(
227  Form("hMCHitMultPerRecoHit_%d", plane),
228  Form("MCHit Multiplicity per RecoHit (Plane %d); Multiplicity; # RecoHits", plane),
229  21,
230  -1.5,
231  19.5));
232  }
233 
235  fs->make<TH1D>("hRecoHitReadTime",
236  "Real Time to Read recob::Hit From Disk; Real Time [ms]; Events",
237  200,
238  0,
239  1000);
240 
242  fs->make<TH1D>("hMCHitReadTime",
243  "Real Time to Read sim::MCHit From Disk; Real Time [ms]; Events",
244  200,
245  0,
246  1000);
247 
248  hAnalysisTime =
249  fs->make<TH1D>("hAnalysisTime",
250  "Analysis Time to Run analyze() Function; Real Time [ms]; Events",
251  200,
252  0,
253  1000);
254 
256  fs->make<TH1D>("hMCHitSearchTime",
257  "Time to Search Multiple sim::MCHit Per RecoHit; RealTime [us]; # RecoHit",
258  200,
259  0,
260  200);
261 
262  hMCHitSearchTimeSum = fs->make<TH1D>(
263  "hMCHitSearchTimeSum",
264  "Time to Search Multiple sim::MCHit for All RecoHit in Event; RealTime [ms]; # Event",
265  200,
266  0,
267  200);
268  }
269 
271 
273  {
274  fAnaWatch.Start();
275 
277  auto const clock_data = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(e);
278 
279  fReadWatch.Start();
280  auto const& mchits_v = *e.getValidHandle<std::vector<sim::MCHitCollection>>(fMCHitModuleName);
281  hMCHitReadTime->Fill(fReadWatch.RealTime() * 1.e3);
282 
283  fReadWatch.Start();
284  auto const& recohits = *e.getValidHandle<std::vector<recob::Hit>>(fRecoHitModuleName);
285  hRecoHitReadTime->Fill(fReadWatch.RealTime() * 1.e3);
286 
287  //
288  // Work on purely MCHit info
289  //
290  std::vector<size_t> mchit_mult(geo->Nplanes(), 0);
291  for (auto const& mchits : mchits_v) {
292 
293  auto plane = geo->ChannelToWire(mchits.Channel()).at(0).Plane;
294 
295  mchit_mult.at(plane) += mchits.size();
296 
297  for (auto const& mchit : mchits)
298 
299  hMCHitQ_v.at(plane)->Fill(mchit.Charge(true));
300  }
301 
302  std::vector<size_t> recohit_mult(geo->Nplanes(), 0);
303 
304  double search_time_sum = 0;
305 
306  // Loop over RecoHit
307  for (size_t hit_index = 0; hit_index < recohits.size(); ++hit_index) {
308 
309  auto const& hit = recohits.at(hit_index);
310 
311  auto const& wire_id = hit.WireID();
312 
313  // Fill purly RecoHit info
314  hRecoHitQ_v.at(wire_id.Plane)->Fill(hit.PeakAmplitude());
315  recohit_mult.at(wire_id.Plane) += 1;
316 
317  // Figure out channel & retrieve MCHitCollection for this channel
318  auto ch = geo->PlaneWireToChannel(wire_id);
319 
320  if (mchits_v.size() <= ch)
321  throw cet::exception(__PRETTY_FUNCTION__)
322  << "Channel " << ch << " not found in MCHit vector!" << std::endl;
323 
324  // Found MCHitCollection
325  auto const& mchits = mchits_v.at(ch);
326 
327  if (ch != mchits.Channel())
328  throw cet::exception(__PRETTY_FUNCTION__)
329  << "MCHit channel & vector index mismatch: " << mchits.Channel() << " != " << ch
330  << std::endl;
331 
332  // Locate corresponding MCHit(s) to this RecoHit
333  sim::MCHit start_time, end_time, peak_time;
334 
335  start_time.SetTime(clock_data.TPCTick2TDC(hit.PeakTimeMinusRMS()), 0);
336  peak_time.SetTime(clock_data.TPCTick2TDC(hit.PeakTime()), 0);
337  end_time.SetTime(clock_data.TPCTick2TDC(hit.PeakTimePlusRMS()), 0);
338 
339  fSearchWatch.Start();
340  auto start_iter = std::lower_bound(mchits.begin(), mchits.end(), start_time);
341  auto end_iter = std::upper_bound(mchits.begin(), mchits.end(), end_time);
342  search_time_sum += fSearchWatch.RealTime() * 1.e3;
343  hMCHitSearchTime->Fill(fSearchWatch.RealTime() * 1.e6);
344 
345  double reco_q = hit.PeakAmplitude();
346  double mc_qsum = 0;
347  double mc_q = 0;
348  double mult = 0;
349  double dt_min = 0;
350  double abs_dt_min = 1e9;
351 
352  // Loop over MCHit(s) that reside in start=>end time of MCHit
353  while (start_iter != end_iter) {
354  mc_qsum += (*start_iter).Charge(true);
355  ++mult;
356  double dt = (*start_iter).PeakTime() - peak_time.PeakTime();
357  double abs_dt = dt;
358  if (abs_dt < 0) abs_dt *= -1;
359 
360  if (abs_dt < abs_dt_min) {
361  abs_dt_min = abs_dt;
362  dt_min = dt;
363  mc_q = (*start_iter).Charge(true);
364  }
365  ++start_iter;
366  }
367 
368  // Fill histograms
369  hMCHitMultPerRecoHit_v.at(wire_id.Plane)->Fill(mult);
370 
371  if (!mc_qsum)
372  hVoidHitQ_v.at(wire_id.Plane)->Fill(reco_q);
373 
374  else {
375 
376  hRecoHitQ_v.at(wire_id.Plane)->Fill(reco_q);
377  hCorrQ_v.at(wire_id.Plane)->Fill(mc_q, reco_q);
378  hCorrQSum_v.at(wire_id.Plane)->Fill(mc_qsum, reco_q);
379  hQSumRatio_v.at(wire_id.Plane)->Fill(reco_q / mc_qsum);
380  hQRatio_v.at(wire_id.Plane)->Fill(reco_q / mc_q);
381  hDT_v.at(wire_id.Plane)->Fill(dt_min);
382  }
383 
384  } // end looping over hits
385 
386  // Fill purely-MCHit and purely-RecoHit multiplicity histograms
387  for (unsigned char plane = 0; plane < geo->Nplanes(); ++plane) {
388  std::cout << mchit_mult.at(plane) << std::endl;
389  hMCHitMult_v.at(plane)->Fill(mchit_mult.at(plane));
390  hRecoHitMult_v.at(plane)->Fill(recohit_mult.at(plane));
391  hCorrMult_v.at(plane)->Fill(mchit_mult.at(plane), recohit_mult.at(plane));
392  }
393 
394  hAnalysisTime->Fill(fAnaWatch.RealTime() * 1.e3);
395  hMCHitSearchTimeSum->Fill(search_time_sum);
396  }
397 }
float PeakTime() const
Getter for start time.
Definition: MCHit.h:70
std::vector< WireID > ChannelToWire(raw::ChannelID_t const channel) const
Returns a list of wires connected to the specified TPC channel.
TH1D * hMCHitSearchTimeSum
Time to search for corresponding MCHit per RecoHit summed over all RecoHits (per event) ...
void SetTime(const float peak, const float width)
Setter function for time.
Definition: MCHit.h:54
void analyze(art::Event const &e) override
Declaration of signal hit object.
MCHitAnaExample(fhicl::ParameterSet const &p)
std::vector< TH1D * > hQRatio_v
Histogram (per plane) for a ratio of RecoHit charge to the closest (in time) RecoHit charge...
std::vector< TH2D * > hCorrMult_v
Histogram (per plane) for # of MCHit vs. # of RecoHit.
EDAnalyzer(fhicl::ParameterSet const &pset)
Definition: EDAnalyzer.cc:6
TH1D * hMCHitReadTime
Time to read MCHit from disk.
std::vector< TH1D * > hDT_v
Histogram (per plane) for time-diff between one Reco hit and the closest (in time) RecoHit...
std::string fMCHitModuleName
Producer module for MCHit.
std::vector< TH1D * > hRecoHitMult_v
Histogram (per plane) for RecoHit multiplicity in each event.
#define DEFINE_ART_MODULE(klass)
Definition: ModuleMacros.h:65
std::vector< TH1D * > hRecoHitQ_v
Histogram (per plane) for RecoHit charge.
std::string fRecoHitModuleName
Producer module for recob::Hit.
T get(std::string const &key) const
Definition: ParameterSet.h:314
std::vector< TH1D * > hQSumRatio_v
std::vector< TH1D * > hMCHitMultPerRecoHit_v
Histogram (per plane) for MCHit multiplicity within start=>end time region of RecoHit.
Detector simulation of raw signals on wires.
std::vector< TH2D * > hCorrQ_v
Histogram (per plane) for charge of one MCHit that is closest (in time) to one RecoHit.
raw::ChannelID_t PlaneWireToChannel(WireID const &wireid) const
Returns the ID of the TPC channel connected to the specified wire.
ValidHandle< PROD > getValidHandle(InputTag const &tag) const
TH1D * hMCHitSearchTime
Time to search for corresponding MCHit per RecoHit.
std::vector< TH1D * > hMCHitQ_v
Histogram (per plane) for all MCHit charge.
TH1D * hRecoHitReadTime
Time to read recob::Hit from disk.
unsigned int Nplanes(TPCID const &tpcid=tpc_zero) const
Returns the total number of planes in the specified TPC.
Definition: GeometryCore.h:977
Float_t e
Definition: plot.C:35
std::vector< TH1D * > hMCHitMult_v
Histogram (per plane) for MCHit multiplicity in each event.
Namespace collecting geometry-related classes utilities.
TH1D * hAnalysisTime
Time to run analyze() function.
art framework interface to geometry description
std::vector< TH2D * > hCorrQSum_v
Histogram (per plane) for charge sum of multiple MCHit found within start=>end time of one RecoHit...
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
std::vector< TH1D * > hVoidHitQ_v
Histogram (per plane) for charge of some RecoHit that has no corresponding MCHit. ...