LArSoft  v08_35_01
Liquid Argon Software toolkit - http://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 
10 #include "art_root_io/TFileService.h"
15 #include "fhiclcpp/ParameterSet.h"
16 
21 #include <TH1D.h>
22 #include <TH2D.h>
23 #include <TStopwatch.h>
24 
25 namespace hit {
26 
27  class MCHitAnaExample;
28 
30  public:
31  explicit MCHitAnaExample(fhicl::ParameterSet const & p);
32  virtual ~MCHitAnaExample();
33 
34  private:
35 
36  void analyze(art::Event const & e) override;
37 
38 
39  //
40  // Module names
41  //
42 
44  std::string fRecoHitModuleName;
45 
47  std::string fMCHitModuleName;
48 
49  //
50  // Stopwatch for timings
51  //
53 
54  //
55  // Histogram for timing record
56  //
57 
60 
63 
66 
69 
72 
73  //
74  // Histograms purely from MCHit
75  //
76 
78  std::vector<TH1D*> hMCHitQ_v;
79 
81  std::vector<TH1D*> hMCHitMult_v;
82 
83  //
84  // Histograms purely from RecoHit
85  //
86 
88  std::vector<TH1D*> hRecoHitQ_v;
89 
91  std::vector<TH1D*> hRecoHitMult_v;
92 
93  //
94  // Histograms from MC+Reco combined/compared information
95  //
96 
98  std::vector<TH2D*> hCorrMult_v;
99 
101  std::vector<TH1D*> hVoidHitQ_v;
102 
104  std::vector<TH2D*> hCorrQ_v;
105 
107  std::vector<TH2D*> hCorrQSum_v;
108 
110  std::vector<TH1D*> hQRatio_v;
111 
116  std::vector<TH1D*> hQSumRatio_v;
117 
119  std::vector<TH1D*> hDT_v;
120 
122  std::vector<TH1D*> hMCHitMultPerRecoHit_v;
123 
124  };
125 
126 
128  :
129  EDAnalyzer(p) // ,
130 
131  {
132  fRecoHitModuleName = p.get<std::string>( "RecoHitModuleName" );
133  fMCHitModuleName = p.get<std::string>( "MCHitModuleName" );
134 
136 
138 
139  for(unsigned char plane=0; plane < geo->Nplanes(); ++plane) {
140 
141  hMCHitQ_v.push_back(fs->make<TH1D>(Form("hMCHitQ_%d",plane),
142  Form("MCHit Charge (Plane %d); Charge [ADC]; # MCHit",plane),
143  250,0,2000));
144 
145  hMCHitMult_v.push_back(fs->make<TH1D>(Form("hMCHitMult_%d",plane),
146  Form("MCHit Multiplicity (Plane %d); Multiplicity; # MCHit",plane),
147  250,0,5000));
148 
149  hRecoHitQ_v.push_back(fs->make<TH1D>(Form("hRecoHitQ_%d",plane),
150  Form("RecoHit Charge (Plane %d); Charge [ADC]; # RecoHits",plane),
151  250,0,2000));
152 
153  hRecoHitMult_v.push_back(fs->make<TH1D>(Form("hRecoHitMult_%d",plane),
154  Form("RecoHit Multiplicity (Plane %d); Multiplicity; # RecoHits",plane),
155  200,0,2000));
156 
157  hCorrMult_v.push_back(fs->make<TH2D>(Form("hCorrMult_%d",plane),
158  Form("# MCHits vs. # RecoHits (Plane %d); # MCHits; # RecoHits",plane),
159  250,0,5000,
160  200,0,2000));
161 
162  hVoidHitQ_v.push_back(fs->make<TH1D>(Form("hVoidHitQ_%d",plane),
163  Form("RecoHit Charge (No Corresponding MCHit)) (Plane %d); Charge [ADC]; # RecoHits",plane),
164  250,0,2000));
165 
166  hCorrQ_v.push_back(fs->make<TH2D>(Form("hCorrQ_%d",plane),
167  Form("RecoHit Q vs. Closest MCHit Q (Plane %d); MCHit Charge [ADC]; RecoHit Charge [ADC]",plane),
168  250,0,2000,
169  250,0,2000));
170 
171  hCorrQSum_v.push_back(fs->make<TH2D>(Form("hCorrQSum_%d",plane),
172  Form("RecoHit Q vs. MCHit QSum Within Start=>End (Plane %d); MCHit Charge [ADC]; RecoHit Charge [ADC]",plane),
173  250,0,2000,
174  250,0,2000));
175 
176  hQRatio_v.push_back(fs->make<TH1D>(Form("hQRatio_%d",plane),
177  Form("Reco/MCHit Charge (Plane %d); Ratio; # RecoHits",plane),
178  200, 0, 2));
179 
180  hQSumRatio_v.push_back(fs->make<TH1D>(Form("hQSumRatio_%d",plane),
181  Form("Reco/MCHit Charge (Plane %d); Ratio; # RecoHits",plane),
182  200, 0, 2));
183 
184  hDT_v.push_back(fs->make<TH1D>(Form("hDT_%d",plane),
185  Form("#DeltaT btw Reco and Closest MCHit (Plane %d); #DeltaT; # RecoHits",plane),
186  190,-9.5,9.5));
187 
188  hMCHitMultPerRecoHit_v.push_back(fs->make<TH1D>(Form("hMCHitMultPerRecoHit_%d",plane),
189  Form("MCHit Multiplicity per RecoHit (Plane %d); Multiplicity; # RecoHits",plane),
190  21,-1.5,19.5));
191  }
192 
193  hRecoHitReadTime = fs->make<TH1D>("hRecoHitReadTime",
194  "Real Time to Read recob::Hit From Disk; Real Time [ms]; Events",
195  200,0,1000);
196 
197  hMCHitReadTime = fs->make<TH1D>("hMCHitReadTime",
198  "Real Time to Read sim::MCHit From Disk; Real Time [ms]; Events",
199  200,0,1000);
200 
201  hAnalysisTime = fs->make<TH1D>("hAnalysisTime",
202  "Analysis Time to Run analyze() Function; Real Time [ms]; Events",
203  200,0,1000);
204 
205  hMCHitSearchTime = fs->make<TH1D>("hMCHitSearchTime",
206  "Time to Search Multiple sim::MCHit Per RecoHit; RealTime [us]; # RecoHit",
207  200,0,200);
208 
209  hMCHitSearchTimeSum = fs->make<TH1D>("hMCHitSearchTimeSum",
210  "Time to Search Multiple sim::MCHit for All RecoHit in Event; RealTime [ms]; # Event",
211  200,0,200);
212 
213  }
214 
216  {}
217 
219  {
220  fAnaWatch.Start();
221 
223  const detinfo::DetectorClocks* ts = lar::providerFrom<detinfo::DetectorClocksService>();
224 
225  fReadWatch.Start();
227  e.getByLabel(fMCHitModuleName,mcHandle);
228  if(!mcHandle.isValid()) throw cet::exception(__PRETTY_FUNCTION__) << "MCHit not found!" <<std::endl;
229  hMCHitReadTime->Fill(fReadWatch.RealTime() * 1.e3);
230  const std::vector<sim::MCHitCollection> &mchits_v = *mcHandle;
231 
232  fReadWatch.Start();
234  e.getByLabel(fRecoHitModuleName,recoHandle);
235  if(!recoHandle.isValid()) throw cet::exception(__PRETTY_FUNCTION__) << "RecoHit not found!" << std::endl;
236  hRecoHitReadTime->Fill(fReadWatch.RealTime() * 1.e3);
237  const std::vector<recob::Hit> &recohits = *recoHandle;
238 
239  //
240  // Work on purely MCHit info
241  //
242  std::vector<size_t> mchit_mult(geo->Nplanes(),0);
243  for(auto const& mchits : mchits_v) {
244 
245  auto plane = geo->ChannelToWire(mchits.Channel()).at(0).Plane;
246 
247  mchit_mult.at(plane) += mchits.size();
248 
249  for(auto const& mchit : mchits)
250 
251  hMCHitQ_v.at(plane)->Fill(mchit.Charge(true));
252 
253  }
254 
255  std::vector<size_t> recohit_mult(geo->Nplanes(),0);
256 
257  double search_time_sum = 0;
258 
259  // Loop over RecoHit
260  for(size_t hit_index=0; hit_index<recohits.size(); ++hit_index) {
261 
262  auto const &hit = recohits.at(hit_index);
263 
264  auto const &wire_id = hit.WireID();
265 
266  // Fill purly RecoHit info
267  hRecoHitQ_v.at(wire_id.Plane)->Fill(hit.PeakAmplitude());
268  recohit_mult.at(wire_id.Plane) += 1;
269 
270  // Figure out channel & retrieve MCHitCollection for this channel
271  auto ch = geo->PlaneWireToChannel(wire_id.Plane,
272  wire_id.Wire,
273  wire_id.TPC,
274  wire_id.Cryostat);
275 
276  if(mchits_v.size() <= ch )
277  throw cet::exception(__PRETTY_FUNCTION__)
278  << "Channel "<<ch<<" not found in MCHit vector!"<<std::endl;
279 
280  // Found MCHitCollection
281  auto const& mchits = mchits_v.at(ch);
282 
283  if( ch != mchits.Channel() )
284  throw cet::exception(__PRETTY_FUNCTION__)
285  << "MCHit channel & vector index mismatch: "
286  << mchits.Channel() << " != " << ch << std::endl;
287 
288  // Locate corresponding MCHit(s) to this RecoHit
289  sim::MCHit start_time, end_time, peak_time;
290 
291  start_time.SetTime ( ts->TPCTick2TDC( hit.PeakTimeMinusRMS() ), 0 );
292  peak_time.SetTime ( ts->TPCTick2TDC( hit.PeakTime() ), 0 );
293  end_time.SetTime ( ts->TPCTick2TDC( hit.PeakTimePlusRMS() ), 0 );
294 
295  fSearchWatch.Start();
296  auto start_iter = std::lower_bound ( mchits.begin(), mchits.end(), start_time );
297  auto end_iter = std::upper_bound ( mchits.begin(), mchits.end(), end_time );
298  search_time_sum += fSearchWatch.RealTime() * 1.e3;
299  hMCHitSearchTime->Fill(fSearchWatch.RealTime() * 1.e6);
300 
301  double reco_q = hit.PeakAmplitude();
302  double mc_qsum = 0;
303  double mc_q = 0;
304  double mult = 0;
305  double dt_min = 0;
306  double abs_dt_min = 1e9;
307 
308  // Loop over MCHit(s) that reside in start=>end time of MCHit
309  while(start_iter != end_iter) {
310  mc_qsum += (*start_iter).Charge(true);
311  ++mult;
312  double dt = (*start_iter).PeakTime() - peak_time.PeakTime();
313  double abs_dt = dt;
314  if(abs_dt<0) abs_dt *= -1;
315 
316  if(abs_dt < abs_dt_min) {
317  abs_dt_min = abs_dt;
318  dt_min = dt;
319  mc_q = (*start_iter).Charge(true);
320  }
321  ++start_iter;
322  }
323 
324  // Fill histograms
325  hMCHitMultPerRecoHit_v.at (wire_id.Plane) -> Fill( mult );
326 
327  if(!mc_qsum) hVoidHitQ_v.at (wire_id.Plane) -> Fill( reco_q );
328 
329  else {
330 
331  hRecoHitQ_v.at (wire_id.Plane) -> Fill ( reco_q );
332  hCorrQ_v.at (wire_id.Plane) -> Fill ( mc_q, reco_q );
333  hCorrQSum_v.at (wire_id.Plane) -> Fill ( mc_qsum, reco_q );
334  hQSumRatio_v.at (wire_id.Plane) -> Fill ( reco_q / mc_qsum );
335  hQRatio_v.at (wire_id.Plane) -> Fill ( reco_q / mc_q );
336  hDT_v.at (wire_id.Plane) -> Fill ( dt_min );
337  }
338 
339  } // end looping over hits
340 
341  // Fill purely-MCHit and purely-RecoHit multiplicity histograms
342  for(unsigned char plane=0; plane<geo->Nplanes(); ++plane) {
343  std::cout<<mchit_mult.at(plane)<<std::endl;
344  hMCHitMult_v.at (plane) -> Fill( mchit_mult.at (plane) );
345  hRecoHitMult_v.at (plane) -> Fill( recohit_mult.at (plane) );
346  hCorrMult_v.at (plane) -> Fill( mchit_mult.at (plane), recohit_mult.at (plane) );
347  }
348 
349  hAnalysisTime->Fill(fAnaWatch.RealTime() * 1.e3);
350  hMCHitSearchTimeSum->Fill(search_time_sum);
351 
352  }
353 }
float PeakTime() const
Getter for start time.
Definition: MCHit.h:75
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:57
void analyze(art::Event const &e) override
Declaration of signal hit object.
TNtupleSim Fill(f1, f2, f3, f4)
MCHitAnaExample(fhicl::ParameterSet const &p)
std::vector< geo::WireID > ChannelToWire(raw::ChannelID_t const channel) const
Returns a list of wires connected to the specified TPC channel.
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.h:27
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...
bool isValid() const
Definition: Handle.h:183
std::string fMCHitModuleName
Producer module for MCHit.
bool getByLabel(std::string const &label, std::string const &instance, Handle< PROD > &result) const
Definition: DataViewImpl.h:437
std::vector< TH1D * > hRecoHitMult_v
Histogram (per plane) for RecoHit multiplicity in each event.
unsigned int Nplanes(unsigned int tpc=0, unsigned int cstat=0) const
Returns the total number of wire planes in the specified TPC.
#define DEFINE_ART_MODULE(klass)
Definition: ModuleMacros.h:56
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:231
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.
Conversion of times between different formats and references.
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.
virtual double TPCTick2TDC(double tick) const =0
Converts a TPC time tick into a electronics time tick.
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.
Float_t e
Definition: plot.C:34
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. ...