Monopole_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // Class: Monopole
3 // Module Type: analyzer
4 // File: Monopole_module.cc
5 //
6 // Generated at Tue Jun 9 10:22:10 2015 by Martin Frank using artmod
7 // from cetpkgsupport v1_08_05.
8 ////////////////////////////////////////////////////////////////////////
9 
11 // #include <canvas/Persistency/Common/FindManyP.h>
13 #include <art/Framework/Core/FindManyP.h>
21 // #include <canvas/Utilities/InputTag.h>
22 #include <art/Utilities/InputTag.h>
23 #include <fhiclcpp/ParameterSet.h>
26 
28 
29 #include "Geometry/Geometry.h"
31 #include "MCCheater/BackTracker.h"
32 #include "Monopole/Constants.h"
33 #include "Monopole/Track3D.h"
34 #include "RawData/RawTrigger.h"
35 #include "RecoBase/CellHit.h"
36 #include "RecoBase/Cluster.h"
37 #include "RecoBase/HoughResult.h"
38 #include "Simulation/Particle.h"
39 // #include "nusimdata/SimulationBase/MCParticle.h"
40 // #include "nusimdata/SimulationBase/MCTruth.h"
42 #include "SimulationBase/MCTruth.h"
43 
44 #include <TH1.h>
45 #include <TTree.h>
46 
47 #include <algorithm>
48 #include <bitset>
49 #include <map>
50 #include <vector>
51 
52 namespace mono
53 {
54  class Monopole;
55 }
56 
58 {
59 public:
60  explicit Monopole(fhicl::ParameterSet const & p);
61 
62  Monopole(Monopole const &) = delete;
63  Monopole(Monopole &&) = delete;
64  Monopole & operator = (Monopole const &) = delete;
65  Monopole & operator = (Monopole &&) = delete;
66 
67  void beginJob() override;
68  void analyze(art::Event const & e) override;
69  void respondToOpenInputFile(art::FileBlock const& fb) override;
70 
71 private:
72 
73  int extent_cell(std::vector<sim::FLSHit> const& hits) const;
74  int extent_plane(std::vector<sim::FLSHit> const& hits) const;
75  TVector3 get_position(sim::FLSHit const& hit) const;
76  unsigned n_mc_hits(mono::Track3D const& track) const;
77  unsigned n_slices(std::vector<rb::Cluster> const& slices) const;
78  double phi(sim::Particle const* particle) const;
79  void split_by_view
80  (std::map<std::string, std::vector<sim::FLSHit> > & hitmap) const;
81  double theta(sim::Particle const* particle) const;
82  geo::View_t view(unsigned plane) const;
83 
84 
85  void tset(std::string const& branch_name, double const& value);
86 
88 
91  TTree *tree_;
92  std::map<std::string, TH1*> h_;
93  std::map<std::string, double> t_;
97 };
98 
99 
100 
101 void mono::Monopole::tset(std::string const& branch_name, double const& value)
102 {
103  auto it = t_.find(branch_name);
104  if (it == t_.end())
106  << __FILE__ << ":" << __LINE__ << "\n"
107  << "Branch " << branch_name << " does not exist!"
108  << std::endl;
109 
110  it->second = value;
111 }
112 
113 
114 
116  EDAnalyzer(p),
117  current_file_name_("not set"),
118  branch_default_value_(-9999),
120 {
121 }
122 
123 
124 
126 {
128 
129  std::cout << "mono::Monopole::respondToOpenInputFile: "
130  << "The current file name is "
132  << std::endl;
133 }
134 
135 
136 
138 {
139  tree_ = tfs_->make<TTree>("Event", "Monopole Event Tree");
140 
141  std::vector<std::string> branch_names =
142  { "run_number", "subrun_number", "event_number", "event_length", "n_hits",
143  "n_clusters", "cluster_n_hits", "cluster_n_hits_x", "cluster_n_hits_y",
144  "reco_n_tracks", "is_mc", "n_slices_monopole", "n_slices_higheremoved",
145  "n_hough_peaks_x_max", "n_hough_peaks_y_max",
146  "slice_sum_adc_max",
147  "gen_dxdz", "gen_dydz", "gen_beta",
148  "mc_dxdz", "mc_dydz", "mc_dxdt", "mc_dydt", "mc_dzdt",
149  "mc_velocity", "mc_beta",
150  "mc_length", "mc_phi", "mc_theta",
151  "mc_center_average_edep", "mc_center_n_hits",
152  "mc_near_end_average_edep", "mc_near_end_n_hits",
153  "mc_far_end_average_edep", "mc_far_end_n_hits",
154  "mc_monopole_hit_detector",
155  "mc_dx", "mc_dy", "mc_dz",
156  "mc_dplane_all", "mc_dplane_x", "mc_dplane_y",
157  "mc_dcell_all", "mc_dcell_x", "mc_dcell_y",
158  "mc_n_hits_all", "mc_n_hits_x", "mc_n_hits_y",
159  "ddt_trigger_decisions",
160  "n_hits_ge_100_ADC",
161  "n_hits_ge_150_ADC",
162  "mc_n_hits_all_without_noise",
163  "mc_n_hits_ge_100_ADC",
164  "mc_n_hits_ge_150_ADC",
165  "mc_n_hits_in_slice4d",
166  "mc_n_hits_in_slice4d_ge_100_ADC",
167  "mc_n_hits_on_cosmic_track",
168  "mc_n_hits_on_cosmic_track_ge_100_ADC",
169  "mc_adc_total",
170  "mc_adc_mean",
171  "mc_adc_variance",
172  "mc_adc_min",
173  "mc_adc_max",
174  "mc_dEdx",
175  "mc_dEdx_Birks",
176  "mc_n_cerenkov",
177  "mc_path_length_sum"
178  };
179 
180  std::vector<std::string> track_names =
181  { "dxdz", "dydz", "dxdt", "dydt", "dzdt", "velocity", "beta",
182  "dzdt_x", "dzdt_y", "r2_xt", "r2_yt", "r2_zt_x", "r2_zt_y",
183  "n_hits", "n_hits_x", "n_hits_y", "n_mc_hits", "sum_adc",
184  "dt", "dx", "dy", "dz", "length", "length_xz", "length_yz",
185  "tstart", "tend",
186  "max_time_gap_xz", "max_time_gap_yz",
187  "most_missing_planes_xz", "most_missing_planes_yz",
188  "dcell_x", "dcell_y", "dplane_x", "dplane_y"
189  };
190 
191  for (auto const& tname : track_names)
192  for (unsigned n_track = 1; n_track <= n_tracks_to_record_; ++n_track)
193  branch_names.
194  push_back("reco_track_" + std::to_string(n_track) + "_" + tname);
195 
196  for (auto const& bname : branch_names)
197  tree_->Branch(bname.c_str(), &t_[bname], (bname + "/D").c_str());
198 
199  tree_->Branch("input_file_name", &current_file_name_);
200 
201 
202  h_["adc_mc"] =
203  tfs_->make<TH1F>("adc_mc", "ADC Distribution - MC", 4096, -0.5, 4095.5);
204  h_["adc_data"] =
205  tfs_->make<TH1F>("adc_data", "ADC Distribution - Data", 4096, -0.5, 4095.5);
206 
207 
208  // 64,000 TDC = 1000 us
209  h_["tdc_mc"] =
210  tfs_->make<TH1F>("tdc_mc", "TDC Distribution - MC", 1000, 0, 64000);
211 
212  h_["dEdx"] =
213  tfs_->make<TH1F>("dEdx", "dE/dx Distribution;dE/dx [MeV/cm]", 200, 0, 40);
214 }
215 
216 
217 
219 {
220  //
221  // Extract Objects from the Event
222  //
223  // Raw Trigger Header
225  e.getByLabel("daq", raw_triggers);
226  if (raw_triggers->size() != 1)
228  << __FILE__ << ":" << __LINE__ << "\n"
229  << "RawTrigger vector has incorrect size: " << raw_triggers->size()
230  << std::endl;
231  unsigned long long
232  event_time(raw_triggers->front().fTriggerTimingMarker_TimeStart);
233  uint32_t event_length(raw_triggers->front().fTriggerRange_TriggerLength);
234 
235  // Offline Cell Hits
237  e.getByLabel("calhit", offline_hits);
238  std::vector<art::Ptr<rb::CellHit> > offline_hitptrs;
239  for (unsigned i = 0; i < offline_hits->size(); ++i)
240  offline_hitptrs.emplace_back(offline_hits, i);
241 
242  /*
243  EXCLUDE EXTRA INFO PARSING FOR NOW
244 
245  uint32_t extra_info(raw_triggers->front().fTriggerHeader_SourceSubID);
246  std::bitset<novaddt::smt::Extra_Info::BITSET_SIZE> extra_bits(extra_info);
247 
248  std::cout << "MF: extra_info = " << extra_info << std::endl;
249  std::cout << "MF: extra_bits = " << extra_bits << std::endl;
250  if (extra_bits.test(novaddt::smt::Extra_Info::INVALID_OLD_DEFAULT_BIT))
251  std::cout << "MF: extra_info is invalid (i.e. set to old default)"
252  << std::endl;
253  else
254  std::cout << "MF: (Check, Timeout, Invalid) = (" << std::boolalpha
255  << extra_bits.test(novaddt::smt::Extra_Info::TIMEOUT_BIT)
256  << ", " << extra_bits.test(novaddt::smt::Extra_Info::CHECK_BIT)
257  << ", " << extra_bits.
258  test(novaddt::smt::Extra_Info::INVALID_OLD_DEFAULT_BIT)
259  << ")" << std::endl;
260  */
261 
262  // DDT Decisions
264  e.getByLabel("slowmonopoletrigger", ddt_decisions);
265 
266  // NOvA Slicer (Slicer4D)
268  e.getByLabel("slicer", nova_slices);
269 
270  // Cosmic Tracks
271  art::Handle<std::vector<rb::Track> > cosmic_tracks;
272  e.getByLabel("cosmictrack", cosmic_tracks);
273 
274  // Monopole Cluster
275  art::Handle<std::vector<rb::Cluster> > mono_clusters;
276  e.getByLabel("monocluster", mono_clusters);
277 
278  // Monopole Slicer
280  e.getByLabel("monoslicer", mono_slices);
281 
282  // High Energy Removed Slices
283  art::Handle<std::vector<rb::Cluster> > higheremoved_slices;
284  e.getByLabel("higheremoval", higheremoved_slices);
285 
286  // Reconstructed Tracks
288  e.getByLabel("monotrack", reco_tracks);
289 
290  // Generator Information
292  e.getByLabel("generator", gen_info);
293 
294 
295  //
296  // Print Hits
297  //
298  // Monopole Track Hits
299  // std::cout << "MF: N_tracks = " << reco_tracks->size() << std::endl;
300  // if (reco_tracks->size() == 1)
301  // std::cout << reco_tracks->front() << std::endl;
302 
303  // // All Hits
304  // art::ServiceHandle<geo::Geometry> geometry;
305  // for (auto const& hit : offline_hitptrs)
306  // {
307  // std::string view;
308  // if (hit->View() == geo::View_t::kX)
309  // view = "XZ";
310  // else if (hit->View() == geo::View_t::kY)
311  // view = "YZ";
312  // else
313  // view = "BAD";
314 
315  // double xyz[3];
316  // geometry->CellInfo(hit->Plane(), hit->Cell(), NULL, xyz);
317  // TVector3 hit_vector(xyz);
318  // double v = geometry->CellTpos(hit->Plane(), hit->Cell());
319  // double z = hit_vector.z();
320  // double t = hit->TNS();
321 
322  // if (z > 1050 && z < 1400 &&
323  // t > -16000 && t < -9000)
324  // std::cout << "AllHits: View " << view
325  // << " Hit: t, v, z = " << t << ", " << v << ", " << z << "\n";
326  // }
327 
328 
329 
330  //
331  // Fill Histograms (and add some hit-based quantities together)
332  //
333  int n_hits_ge_100_ADC = 0;
334  int n_hits_ge_150_ADC = 0;
335 
336  double mc_n_hits_all_without_noise = 0.0;
337  int mc_n_hits_ge_100_ADC = 0;
338  int mc_n_hits_ge_150_ADC = 0;
339 
340  double mc_adc_total = 0.0;
341  double mc_adc_mean = 0.0;
342  double mc_adc_variance = 0.0;
343  double mc_adc_max = -9999;
344  double mc_adc_min = -9999;
345  if (!offline_hits.failedToGet())
346  {
347  std::vector<int32_t> mc_tdc_values;
348  std::vector<double> mc_adc_values;
349  for (auto const& hit : offline_hitptrs)
350  {
351  if (hit->ADC() >= 100)
352  ++n_hits_ge_100_ADC;
353 
354  if (hit->ADC() >= 150)
355  ++n_hits_ge_150_ADC;
356 
357  if (hit->IsMC())
358  {
359  if (!bt_->IsNoise(hit))
360  {
361  ++mc_n_hits_all_without_noise;
362 
363  h_.at("adc_mc")->Fill(hit->ADC());
364  mc_adc_values.push_back(hit->ADC());
365  mc_tdc_values.push_back(hit->TDC());
366 
367  if (hit->ADC() >= 100)
368  ++mc_n_hits_ge_100_ADC;
369 
370  if (hit->ADC() >= 150)
371  ++mc_n_hits_ge_150_ADC;
372  }
373  } else {
374  h_.at("adc_data")->Fill(hit->ADC());
375  }
376  }
377 
378  if (!mc_tdc_values.empty())
379  {
380  auto tdc_min =
381  *std::min_element(mc_tdc_values.cbegin(), mc_tdc_values.cend());
382  for (auto const& tdc : mc_tdc_values)
383  h_.at("tdc_mc")->Fill(tdc - tdc_min);
384  }
385 
386  if (mc_n_hits_all_without_noise > 1)
387  {
388  for (auto const& adc : mc_adc_values)
389  mc_adc_total += adc;
390 
391  mc_adc_mean = mc_adc_total / mc_n_hits_all_without_noise;
392 
393  for (auto const& adc : mc_adc_values)
394  mc_adc_variance += std::pow(adc - mc_adc_mean, 2);
395 
396  mc_adc_variance = std::sqrt(mc_adc_variance /
397  (mc_n_hits_all_without_noise - 1));
398 
399  mc_adc_min =
400  *std::min_element(mc_adc_values.cbegin(), mc_adc_values.cend());
401  mc_adc_max =
402  *std::max_element(mc_adc_values.cbegin(), mc_adc_values.cend());
403  }
404  }
405 
406  // Slicer4D Hit Counts
407  int mc_n_hits_in_slice4d = 0;
408  int mc_n_hits_in_slice4d_ge_100_ADC = 0;
409  if (!nova_slices.failedToGet())
410  {
411  for (auto const& nova_slice : *nova_slices)
412  {
413  if (!nova_slice.IsNoise())
414  {
415  for (auto const& nova_slice_hit : nova_slice.AllCells())
416  {
417  if (nova_slice_hit->IsMC() &&
418  !bt_->IsNoise(nova_slice_hit))
419  {
420  ++mc_n_hits_in_slice4d;
421 
422  if (nova_slice_hit->ADC() >= 100)
423  ++mc_n_hits_in_slice4d_ge_100_ADC;
424  }
425  }
426  }
427  }
428  }
429 
430  // Cosmic Track Hit Counts
431  int mc_n_hits_on_cosmic_track = 0;
432  int mc_n_hits_on_cosmic_track_ge_100_ADC = 0;
433  if (!cosmic_tracks.failedToGet())
434  {
435  for (auto const& cosmic_track : *cosmic_tracks)
436  {
437  for (auto const& cosmic_track_hit : cosmic_track.AllCells())
438  {
439  if (cosmic_track_hit->IsMC() &&
440  !bt_->IsNoise(cosmic_track_hit))
441  {
442  ++mc_n_hits_on_cosmic_track;
443 
444  if (cosmic_track_hit->ADC() >= 100)
445  ++mc_n_hits_on_cosmic_track_ge_100_ADC;
446  }
447  }
448  }
449  }
450 
451 
452 
453  //
454  // Fill Tree
455  //
456  for (auto & branch : t_)
457  branch.second = branch_default_value_;
458 
459  // Event Information
460  tset("run_number", e.run());
461  tset("subrun_number", e.subRun());
462  tset("event_number", e.event());
463  tset("event_length", event_length);
464 
465  // DDT Information
466  if (!ddt_decisions.failedToGet())
467  tset("ddt_trigger_decisions", ddt_decisions->size());
468 
469  // Hit-Level Information (from histogram filling above)
470  tset("n_hits_ge_100_ADC", n_hits_ge_100_ADC);
471  tset("n_hits_ge_150_ADC", n_hits_ge_150_ADC);
472 
473  tset("mc_n_hits_all_without_noise", mc_n_hits_all_without_noise);
474  tset("mc_n_hits_ge_100_ADC", mc_n_hits_ge_100_ADC);
475  tset("mc_n_hits_ge_150_ADC", mc_n_hits_ge_150_ADC);
476 
477  tset("mc_adc_total", mc_adc_total);
478  tset("mc_adc_mean", mc_adc_mean);
479  tset("mc_adc_variance", mc_adc_variance);
480  tset("mc_adc_min", mc_adc_min);
481  tset("mc_adc_max", mc_adc_max);
482 
483  // NOvA Slice (Slicer4D) Information
484  tset("mc_n_hits_in_slice4d", mc_n_hits_in_slice4d);
485  tset("mc_n_hits_in_slice4d_ge_100_ADC", mc_n_hits_in_slice4d_ge_100_ADC);
486 
487  // Cosmic Track Information
488  tset("mc_n_hits_on_cosmic_track", mc_n_hits_on_cosmic_track);
489  tset("mc_n_hits_on_cosmic_track_ge_100_ADC",
490  mc_n_hits_on_cosmic_track_ge_100_ADC);
491 
492  // Reconstruction Information
493  if (!offline_hits.failedToGet())
494  tset("n_hits", offline_hits->size());
495 
496  if (!mono_clusters.failedToGet())
497  {
498  tset("n_clusters", mono_clusters->size());
499 
500  if (mono_clusters->size() == 1)
501  {
502  auto cluster = mono_clusters->front();
503  tset("cluster_n_hits", cluster.NCell());
504  tset("cluster_n_hits_x", cluster.NXCell());
505  tset("cluster_n_hits_y", cluster.NYCell());
506  }
507  }
508 
509  if (!mono_slices.failedToGet())
510  {
511  tset("n_slices_monopole", n_slices(*mono_slices));
512 
513  auto max_it = std::max_element
514  (mono_slices->begin(), mono_slices->end(),
515  [](auto const lhs, auto const rhs)
516  { return lhs.TotalADC() < rhs.TotalADC(); });
517  tset("slice_sum_adc_max", max_it->TotalADC());
518  }
519 
520  if (!higheremoved_slices.failedToGet())
521  {
522  tset("n_slices_higheremoved", n_slices(*higheremoved_slices));
523 
524  // Retrieve Hough information from the event:
526  find_hough_results(higheremoved_slices, e, "hough");
527 
528  // Record the maximum number of Hough peaks:
529  std::map<geo::View_t, unsigned> n_peaks_max =
530  { {geo::View_t::kX, 0},
531  {geo::View_t::kY, 0} };
532 
533  for (size_t n_slice = 0; n_slice != higheremoved_slices->size(); ++n_slice)
534  {
535  if (!higheremoved_slices->at(n_slice).IsNoise())
536  {
537  std::vector<art::Ptr<rb::HoughResult> > hough_results;
538  find_hough_results.get(n_slice, hough_results);
539 
540  for (auto const& hough_result : hough_results)
541  {
542  auto view = hough_result->fView;
543  unsigned n_peak = hough_result->fPeak.size();
544  n_peaks_max.at(view) = std::max(n_peak, n_peaks_max.at(view));
545  }
546  }
547  }
548  tset("n_hough_peaks_x_max", n_peaks_max.at(geo::View_t::kX));
549  tset("n_hough_peaks_y_max", n_peaks_max.at(geo::View_t::kY));
550  }
551 
552  if (!reco_tracks.failedToGet())
553  {
554  tset("reco_n_tracks", reco_tracks->size());
555 
556  unsigned n_tracks =
557  std::min(n_tracks_to_record_, unsigned(reco_tracks->size()));
558  for (unsigned n_track = 1; n_track <= n_tracks; ++n_track)
559  {
560  auto track = reco_tracks->at(n_track - 1);
561 
562  std::string name = "reco_track_" + std::to_string(n_track) + "_";
563  tset(name + "dxdz", track.dxdz());
564  tset(name + "dydz", track.dydz());
565 
566  tset(name + "dxdt", track.velocity().x());
567  tset(name + "dydt", track.velocity().y());
568  tset(name + "dzdt", track.velocity().z());
569  tset(name + "velocity", track.velocity().Mag());
570  tset(name + "beta", track.beta());
571 
572  tset(name + "dzdt_x", track.dzdt(geo::View_t::kX));
573  tset(name + "dzdt_y", track.dzdt(geo::View_t::kY));
574 
575  // squared linear correlation factors:
576  tset(name + "r2_xt", track.r2_xt());
577  tset(name + "r2_yt", track.r2_yt());
578  tset(name + "r2_zt_x", track.r2_zt(geo::View_t::kX));
579  tset(name + "r2_zt_y", track.r2_zt(geo::View_t::kY));
580 
581  tset(name + "n_hits", track.cluster().NCell());
582  tset(name + "n_hits_x", track.cluster().NXCell());
583  tset(name + "n_hits_y", track.cluster().NYCell());
584  tset(name + "n_mc_hits", n_mc_hits(track));
585 
586  tset(name + "sum_adc", track.cluster().TotalADC());
587 
588  tset(name + "tstart", track.cluster().MinTNS());
589  tset(name + "tend", track.cluster().MaxTNS());
590 
591  tset(name + "dcell_x", track.cluster().ExtentCell(geo::View_t::kX));
592  tset(name + "dcell_y", track.cluster().ExtentCell(geo::View_t::kY));
593 
594  tset(name + "dplane_x", track.cluster().ExtentPlane(geo::View_t::kX));
595  tset(name + "dplane_y", track.cluster().ExtentPlane(geo::View_t::kY));
596 
597  TVector3 ext = track.cluster().ExtentXYZ();
598  tset(name + "dx", ext.x());
599  tset(name + "dy", ext.y());
600  tset(name + "dz", ext.z());
601  tset(name + "length", ext.Mag());
602  tset(name + "length_xz", util::pythag(ext.x(), ext.z()));
603  tset(name + "length_yz", util::pythag(ext.y(), ext.z()));
604  tset(name + "dt", track.cluster().ExtentTNS());
605 
606  tset(name + "max_time_gap_xz", track.max_time_gap_xz());
607  tset(name + "max_time_gap_yz", track.max_time_gap_yz());
608 
609  tset(name + "most_missing_planes_xz",
610  track.cluster().MostMissingPlanes(geo::View_t::kX));
611  tset(name + "most_missing_planes_yz",
612  track.cluster().MostMissingPlanes(geo::View_t::kY));
613  }
614  }
615 
616  if (bt_->HaveTruthInfo())
617  {
618  tset("is_mc", true);
619  tset("mc_monopole_hit_detector", false);
620 
622 
623  // Generator Information
624  /*
625  While working with the generator information, I found that its
626  information is incorrect (e.g. m = -2e9 GeV/c^2 for monopole).
627  The information extracted from the particle navigator, however,
628  is correct.
629 
630  TO DO:
631  1. since the generator mass is incorrect, implement it manually:
632  const double MONOPOLE_MASS = 1e16; // GeV / c^2
633 
634  2. use the hard-coded mass to calculate all of the generator quantities
635 
636  3. add all of the MC quantities at the generator level
637 
638  4. add the mass, px, py, pz, and p to MC and Gen level information
639  (makes it easy to calculate all quantities with the tree)
640  */
641  if (gen_info->size() == 1)
642  {
643  simb::MCParticle const& particle = gen_info->front().GetParticle(0);
644  unsigned const trajectory = 0;
645 
646  auto m = particle.Mass();
647  auto p = particle.Momentum(trajectory);
648 
649  tset("gen_dxdz", p.Px() / p.Pz());
650  tset("gen_dydz", p.Py() / p.Pz());
651  tset("gen_beta", p.P() / m);
652  } else {
654  << __FILE__ << ":" << __LINE__ << "\n"
655  << "The generator object contains " << gen_info->size()
656  << " elements, there should only be 1." << std::endl;
657  }
658 
659  // MC Information
661  std::cout << "MF: Number of MC Particles in Navigator = "
662  << pn.size() << std::endl;
663  if (!pn.empty())
664  {
665  tset("mc_monopole_hit_detector", true);
666 
667  sim::Particle const* particle = pn.Primary(0);
668  unsigned const trajectory = 0;
669 
670  auto m = particle->Mass();
671  auto p = particle->Momentum(trajectory);
672 
673  tset("mc_phi", phi(particle));
674  tset("mc_theta", theta(particle));
675 
676  tset("mc_dxdz", p.Px() / p.Pz());
677  tset("mc_dydz", p.Py() / p.Pz());
678 
679  tset("mc_dxdt", c * p.Px() / m);
680  tset("mc_dydt", c * p.Py() / m);
681  tset("mc_dzdt", c * p.Pz() / m);
682  tset("mc_velocity", c * p.P() / m);
683  tset("mc_beta", p.P() / m);
684 
685  std::map<std::string, std::vector<sim::FLSHit> > fls_hits;
686  fls_hits["all"] = bt_->ParticleToFLSHit(particle->TrackId());
687  split_by_view(fls_hits);
688 
689  // dE/dx
690  double Edep_sum = 0.0;
691  double Edep_Birks_sum = 0.0;
692  double path_length_sum = 0.0;
693  double n_cerenkov_sum = 0.0;
694  for (auto hit : fls_hits["all"])
695  {
696  double path_length = hit.GetTotalPathLength();
697  if (path_length > 5.0)
698  {
699  Edep_sum += hit.GetEdep();
700  Edep_Birks_sum += hit.GetEdepBirks();
701  path_length_sum += path_length;
702 
703  h_.at("dEdx")->Fill(hit.GetEdep() / path_length * 1e3); // MeV/cm
704  }
705  }
706 
707  tset("mc_path_length_sum", path_length_sum); // cm
708  if (path_length_sum > 0)
709  {
710  tset("mc_dEdx", Edep_sum / path_length_sum); // GeV/cm
711  tset("mc_dEdx_Birks", Edep_Birks_sum / path_length_sum); // GeV/cm
712  tset("mc_n_cerenkov", n_cerenkov_sum / path_length_sum);
713  }
714 
715  for (auto const& hits : fls_hits)
716  {
717  std::string suffix = hits.first;
718 
719  tset("mc_n_hits_" + suffix, hits.second.size());
720  tset("mc_dplane_" + suffix, extent_plane(hits.second));
721  tset("mc_dcell_" + suffix, extent_cell(hits.second));
722  }
723 
724  std::sort(fls_hits.at("all").begin(), fls_hits.at("all").end(),
725  [](sim::FLSHit const& lhs, sim::FLSHit const& rhs)
726  { return lhs.GetEntryT() < rhs.GetEntryT(); });
727 
728  auto entry_hit = fls_hits.at("all").front();
729  auto exit_hit = fls_hits.at("all").back();
730  TVector3 entry_position = get_position(entry_hit);
731  TVector3 exit_position = get_position(exit_hit);
732  TVector3 track = exit_position - entry_position;
733  tset("mc_dx", track.X());
734  tset("mc_dy", track.Y());
735  tset("mc_dz", track.Z());
736  tset("mc_length", track.Mag());
737  }
738 
739  // Energy Deposition
740  std::map<std::string, std::vector<double> > edeps;
741  for (auto const& hit : offline_hitptrs)
742  {
743  auto fls_hits = bt_->HitToFLSHit(hit);
744  for (auto const& fls_hit : fls_hits)
745  {
746  double edep =
747  static_cast<double>(hit->ADC()) / fls_hit.GetTotalPathLength();
748 
749  double z = (fls_hit.GetEntryZ() + fls_hit.GetExitZ()) / 2.0;
750  if (-25.0 < z && z < 25.0)
751  edeps["center"].push_back(edep);
752  else if (-650.0 < z && z < -600.0)
753  edeps["far_end"].push_back(edep);
754  else if (600.0 < z && z < 650)
755  edeps["near_end"].push_back(edep);
756  }
757  }
758 
759  for (auto const& position : edeps)
760  {
761  double sum = 0.0;
762  for (auto const& value : position.second)
763  sum += value;
764 
765  double average_edep = sum / static_cast<double>(position.second.size());
766  tset("mc_" + position.first + "_average_edep", average_edep);
767  tset("mc_" + position.first + "_n_hits", position.second.size());
768  }
769  } else {
770  tset("is_mc", false);
771  }
772 
773  tree_->Fill();
774 }
775 
776 
777 
779 (std::vector<sim::FLSHit> const& hits) const
780 {
781  if (hits.empty())
782  return branch_default_value_;
783 
784  auto it = std::minmax_element
785  (hits.cbegin(), hits.cend(),
786  [](auto const lhs, auto const rhs)
787  { return lhs.GetCellID() < rhs.GetCellID(); });
788 
789  // we add 1 to be consistent with rb::Cluster::ExtentCell()
790  return std::abs(it.first->GetCellID() - it.second->GetCellID()) + 1;
791 }
792 
793 
794 
796 (std::vector<sim::FLSHit> const& hits) const
797 {
798  if (hits.empty())
799  return branch_default_value_;
800 
801  auto it = std::minmax_element
802  (hits.cbegin(), hits.cend(),
803  [](auto const lhs, auto const rhs)
804  { return lhs.GetPlaneID() < rhs.GetPlaneID(); });
805 
806  // we add 1 to be consistent with rb::Cluster::ExtentPlane()
807  return std::abs(it.first->GetPlaneID() - it.second->GetPlaneID()) + 1;
808 }
809 
810 
811 
813 {
814  TVector3 result;
815  geo_->Plane(hit.GetPlaneID())->Cell(hit.GetCellID())->
816  GetCenter(result, hit.GetZAverage());
817 
818  return result;
819 }
820 
821 
822 
824 {
825  unsigned result = 0;
826 
827  if (!bt_->HaveTruthInfo())
828  return result;
829 
830  for (auto const& hit : track.cluster().AllCells())
831  if (hit->IsMC())
832  ++result;
833 
834  return result;
835 }
836 
837 
838 
839 unsigned mono::Monopole::n_slices(std::vector<rb::Cluster> const& slices) const
840 {
841  unsigned result = std::count_if
842  (slices.begin(), slices.end(),
843  [](rb::Cluster const& slice){ return !slice.IsNoise(); });
844 
845  return result;
846 }
847 
848 
849 
850 double mono::Monopole::phi(sim::Particle const* particle) const
851 {
852  double Px(particle->Px()), Py(particle->Py());
853 
854  return std::atan2(Py, Px);
855 }
856 
857 
858 
860 (std::map<std::string, std::vector<sim::FLSHit> > & hitmap) const
861 {
862  hitmap["x"].clear();
863  hitmap["y"].clear();
864 
865  for (auto const & hit : hitmap.at("all"))
866  {
867  if (view(hit.GetPlaneID()) == geo::View_t::kX)
868  hitmap.at("x").push_back(hit);
869  else
870  hitmap.at("y").push_back(hit);
871  }
872 }
873 
874 
875 
876 double mono::Monopole::theta(sim::Particle const* particle) const
877 {
878  double Px(particle->Px()), Py(particle->Py()), Pz(particle->Pz());
879  double Pt = std::sqrt(Px * Px + Py * Py);
880 
881  return std::atan2(Pt, Pz);
882 }
883 
884 
885 
887 {
888  // x planes are odd, y planes are even
889  if (plane % 2 == 1)
890  return geo::View_t::kX;
891 
892  return geo::View_t::kY;
893 }
894 
895 
896 
T max(const caf::Proxy< T > &a, T b)
const XML_Char * name
Definition: expat.h:151
SubRunNumber_t subRun() const
Definition: Event.h:72
back track the reconstruction to the simulation
void tset(std::string const &branch_name, double const &value)
std::string current_file_name_
double Py(const int i=0) const
Definition: MCParticle.h:230
set< int >::iterator it
int GetPlaneID() const
Plane ID.
Definition: FLSHit.h:37
Definition: Cluster.h:13
const sim::ParticleNavigator & ParticleNavigator() const
Get a reference to the ParticleNavigator.
Definition: BackTracker.h:744
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
size_type size() const
art::ServiceHandle< art::TFileService > tfs_
int GetCellID() const
Cell ID.
Definition: FLSHit.h:39
const char * p
Definition: xmltok.h:285
T sqrt(T number)
Definition: d0nt_math.hpp:156
double theta(sim::Particle const *particle) const
std::vector< sim::FLSHit > ParticleToFLSHit(const int &trackID) const
All the FLSHits that were created by the track id trackID, sorted from most to least light...
double Mass() const
Definition: MCParticle.h:238
constexpr T pow(T x)
Definition: pow.h:75
TVector3 get_position(sim::FLSHit const &hit) const
double Px(const int i=0) const
Definition: MCParticle.h:229
unsigned n_slices(std::vector< rb::Cluster > const &slices) const
std::string const & fileName() const
Definition: FileBlock.h:38
Definition: event.h:19
A collection of associated CellHits.
Definition: Cluster.h:47
A single unit of energy deposition in the liquid scintillator.
Definition: FLSHit.h:19
const double SPEED_OF_LIGHT
Definition: Constants.h:8
const PlaneGeo * Plane(unsigned int i) const
DEFINE_ART_MODULE(TestTMapFile)
Particle class.
float abs(float number)
Definition: d0nt_math.hpp:39
const Var kY([](const caf::SRProxy *sr){float tmp=0.f;if(sr->mc.nu.empty()) return tmp;tmp=sr->mc.nu[0].y;return tmp;})
art::PtrVector< rb::CellHit > AllCells() const
Get all cells from both views.
Definition: Cluster.cxx:180
bool IsNoise(const art::Ptr< rb::CellHit > &hit) const
Is this hit not associated with any particles?
int TrackId() const
Definition: MCParticle.h:209
int extent_plane(std::vector< sim::FLSHit > const &hits) const
double branch_default_value_
void hits()
Definition: readHits.C:15
const sim::Particle * Primary(const int) const
const XML_Char int const XML_Char * value
Definition: expat.h:331
void split_by_view(std::map< std::string, std::vector< sim::FLSHit > > &hitmap) const
std::vector< sim::FLSHit > HitToFLSHit(const rb::CellHit &hit) const
All the FLSHits that contributed to this hit, sorted from most to least light.
base_types push_back(int_type())
rb::Cluster cluster() const
Definition: Track3D.cxx:185
Monopole & operator=(Monopole const &)=delete
double phi(sim::Particle const *particle) const
void analyze(art::Event const &e) override
void beginJob() override
EventNumber_t event() const
Definition: Event.h:67
EDAnalyzer(Table< Config > const &config)
Definition: EDAnalyzer.h:100
Double_t edep
Definition: macro.C:13
z
Definition: test.py:28
unsigned n_mc_hits(mono::Track3D const &track) const
OStream cout
Definition: OStream.cxx:6
void respondToOpenInputFile(art::FileBlock const &fb) override
float GetEntryT() const
Definition: FLSHit.h:51
cet::coded_exception< errors::ErrorCodes, ExceptionDetail::translate > Exception
Definition: Exception.h:66
std::map< std::string, double > t_
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
T * make(ARGS...args) const
unsigned n_tracks_to_record_
Data resulting from a Hough transform on the cell hit positions.
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
std::map< std::string, TH1 * > h_
double pythag(double x, double y)
2D Euclidean distance
Definition: MathUtil.h:29
Definition: event.h:1
const TLorentzVector & Momentum(const int i=0) const
Definition: MCParticle.h:219
bool HaveTruthInfo() const
Is this a file with truth info in? (Is BackTracker going to be any use to you?)
Definition: BackTracker.h:133
Definition: geo.h:1
double Pz(const int i=0) const
Definition: MCParticle.h:231
Monopole(fhicl::ParameterSet const &p)
size_type get(size_type i, reference item, data_reference data) const
Definition: FindManyP.h:469
geo::View_t view(unsigned plane) const
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
T min(const caf::Proxy< T > &a, T b)
float GetZAverage(const int step) const
Get Z-average for the step. This is in local coordinates.
Definition: FLSHit.h:80
Double_t sum
Definition: plot.C:31
Float_t e
Definition: plot.C:35
RunNumber_t run() const
Definition: Event.h:77
art::ServiceHandle< geo::Geometry > geo_
art::ServiceHandle< cheat::BackTracker > bt_
Encapsulate the geometry of one entire detector (near, far, ndos)
bool failedToGet() const
Definition: Handle.h:196
T atan2(T number)
Definition: d0nt_math.hpp:72
int extent_cell(std::vector< sim::FLSHit > const &hits) const