15 #include "art_root_io/TFileService.h" 24 #include <TStopwatch.h> 28 class MCHitAnaExample;
134 for (
unsigned char plane = 0; plane < geo->
Nplanes(); ++plane) {
137 fs->make<TH1D>(Form(
"hMCHitQ_%d", plane),
138 Form(
"MCHit Charge (Plane %d); Charge [ADC]; # MCHit", plane),
144 fs->make<TH1D>(Form(
"hMCHitMult_%d", plane),
145 Form(
"MCHit Multiplicity (Plane %d); Multiplicity; # MCHit", plane),
151 fs->make<TH1D>(Form(
"hRecoHitQ_%d", plane),
152 Form(
"RecoHit Charge (Plane %d); Charge [ADC]; # RecoHits", plane),
158 fs->make<TH1D>(Form(
"hRecoHitMult_%d", plane),
159 Form(
"RecoHit Multiplicity (Plane %d); Multiplicity; # RecoHits", plane),
165 fs->make<TH2D>(Form(
"hCorrMult_%d", plane),
166 Form(
"# MCHits vs. # RecoHits (Plane %d); # MCHits; # RecoHits", plane),
175 Form(
"hVoidHitQ_%d", plane),
176 Form(
"RecoHit Charge (No Corresponding MCHit)) (Plane %d); Charge [ADC]; # RecoHits",
183 Form(
"hCorrQ_%d", plane),
184 Form(
"RecoHit Q vs. Closest MCHit Q (Plane %d); MCHit Charge [ADC]; RecoHit Charge [ADC]",
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]",
206 fs->make<TH1D>(Form(
"hQRatio_%d", plane),
207 Form(
"Reco/MCHit Charge (Plane %d); Ratio; # RecoHits", plane),
213 fs->make<TH1D>(Form(
"hQSumRatio_%d", plane),
214 Form(
"Reco/MCHit Charge (Plane %d); Ratio; # RecoHits", plane),
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),
227 Form(
"hMCHitMultPerRecoHit_%d", plane),
228 Form(
"MCHit Multiplicity per RecoHit (Plane %d); Multiplicity; # RecoHits", plane),
235 fs->make<TH1D>(
"hRecoHitReadTime",
236 "Real Time to Read recob::Hit From Disk; Real Time [ms]; Events",
242 fs->make<TH1D>(
"hMCHitReadTime",
243 "Real Time to Read sim::MCHit From Disk; Real Time [ms]; Events",
249 fs->make<TH1D>(
"hAnalysisTime",
250 "Analysis Time to Run analyze() Function; Real Time [ms]; Events",
256 fs->make<TH1D>(
"hMCHitSearchTime",
257 "Time to Search Multiple sim::MCHit Per RecoHit; RealTime [us]; # RecoHit",
263 "hMCHitSearchTimeSum",
264 "Time to Search Multiple sim::MCHit for All RecoHit in Event; RealTime [ms]; # Event",
290 std::vector<size_t> mchit_mult(geo->
Nplanes(), 0);
291 for (
auto const& mchits : mchits_v) {
293 auto plane = geo->
ChannelToWire(mchits.Channel()).at(0).Plane;
295 mchit_mult.at(plane) += mchits.size();
297 for (
auto const& mchit : mchits)
299 hMCHitQ_v.at(plane)->Fill(mchit.Charge(
true));
302 std::vector<size_t> recohit_mult(geo->
Nplanes(), 0);
304 double search_time_sum = 0;
307 for (
size_t hit_index = 0; hit_index < recohits.size(); ++hit_index) {
309 auto const&
hit = recohits.at(hit_index);
311 auto const& wire_id =
hit.WireID();
315 recohit_mult.at(wire_id.Plane) += 1;
320 if (mchits_v.size() <= ch)
322 <<
"Channel " << ch <<
" not found in MCHit vector!" << std::endl;
325 auto const& mchits = mchits_v.at(ch);
327 if (ch != mchits.Channel())
329 <<
"MCHit channel & vector index mismatch: " << mchits.Channel() <<
" != " << ch
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);
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);
345 double reco_q =
hit.PeakAmplitude();
350 double abs_dt_min = 1e9;
353 while (start_iter != end_iter) {
354 mc_qsum += (*start_iter).Charge(
true);
356 double dt = (*start_iter).PeakTime() - peak_time.
PeakTime();
358 if (abs_dt < 0) abs_dt *= -1;
360 if (abs_dt < abs_dt_min) {
363 mc_q = (*start_iter).Charge(
true);
377 hCorrQ_v.at(wire_id.Plane)->Fill(mc_q, reco_q);
378 hCorrQSum_v.at(wire_id.Plane)->Fill(mc_qsum, reco_q);
380 hQRatio_v.at(wire_id.Plane)->Fill(reco_q / mc_q);
381 hDT_v.at(wire_id.Plane)->Fill(dt_min);
387 for (
unsigned char plane = 0; plane < geo->
Nplanes(); ++plane) {
388 std::cout << mchit_mult.at(plane) << std::endl;
391 hCorrMult_v.at(plane)->Fill(mchit_mult.at(plane), recohit_mult.at(plane));
float PeakTime() const
Getter for start time.
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.
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)
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)
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
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.
std::vector< TH1D * > hMCHitMult_v
Histogram (per plane) for MCHit multiplicity in each event.
Namespace collecting geometry-related classes utilities.
virtual ~MCHitAnaExample()
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
std::vector< TH1D * > hVoidHitQ_v
Histogram (per plane) for charge of some RecoHit that has no corresponding MCHit. ...