13 #include <art/Framework/Core/FindManyP.h> 22 #include <art/Utilities/InputTag.h> 74 int extent_plane(std::vector<sim::FLSHit>
const& hits)
const;
77 unsigned n_slices(std::vector<rb::Cluster>
const& slices)
const;
80 (std::map<
std::string, std::vector<sim::FLSHit> > & hitmap)
const;
92 std::map<std::string, TH1*>
h_;
93 std::map<std::string, double>
t_;
103 auto it =
t_.find(branch_name);
106 << __FILE__ <<
":" << __LINE__ <<
"\n" 107 <<
"Branch " << branch_name <<
" does not exist!" 129 std::cout <<
"mono::Monopole::respondToOpenInputFile: " 130 <<
"The current file name is " 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",
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",
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",
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",
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" 191 for (
auto const&
tname : track_names)
196 for (
auto const& bname : branch_names)
197 tree_->Branch(bname.c_str(), &
t_[bname], (bname +
"/D").c_str());
203 tfs_->
make<TH1F>(
"adc_mc",
"ADC Distribution - MC", 4096, -0.5, 4095.5);
205 tfs_->
make<TH1F>(
"adc_data",
"ADC Distribution - Data", 4096, -0.5, 4095.5);
210 tfs_->
make<TH1F>(
"tdc_mc",
"TDC Distribution - MC", 1000, 0, 64000);
213 tfs_->
make<TH1F>(
"dEdx",
"dE/dx Distribution;dE/dx [MeV/cm]", 200, 0, 40);
226 if (raw_triggers->size() != 1)
228 << __FILE__ <<
":" << __LINE__ <<
"\n" 229 <<
"RawTrigger vector has incorrect size: " << raw_triggers->size()
232 event_time(raw_triggers->front().fTriggerTimingMarker_TimeStart);
233 uint32_t event_length(raw_triggers->front().fTriggerRange_TriggerLength);
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);
264 e.
getByLabel(
"slowmonopoletrigger", ddt_decisions);
284 e.
getByLabel(
"higheremoval", higheremoved_slices);
333 int n_hits_ge_100_ADC = 0;
334 int n_hits_ge_150_ADC = 0;
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;
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())
347 std::vector<int32_t> mc_tdc_values;
348 std::vector<double> mc_adc_values;
349 for (
auto const&
hit : offline_hitptrs)
351 if (
hit->ADC() >= 100)
354 if (
hit->ADC() >= 150)
361 ++mc_n_hits_all_without_noise;
363 h_.at(
"adc_mc")->Fill(
hit->ADC());
364 mc_adc_values.push_back(
hit->ADC());
365 mc_tdc_values.push_back(
hit->TDC());
367 if (
hit->ADC() >= 100)
368 ++mc_n_hits_ge_100_ADC;
370 if (
hit->ADC() >= 150)
371 ++mc_n_hits_ge_150_ADC;
374 h_.at(
"adc_data")->Fill(
hit->ADC());
378 if (!mc_tdc_values.empty())
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);
386 if (mc_n_hits_all_without_noise > 1)
388 for (
auto const& adc : mc_adc_values)
391 mc_adc_mean = mc_adc_total / mc_n_hits_all_without_noise;
393 for (
auto const& adc : mc_adc_values)
394 mc_adc_variance +=
std::pow(adc - mc_adc_mean, 2);
396 mc_adc_variance =
std::sqrt(mc_adc_variance /
397 (mc_n_hits_all_without_noise - 1));
400 *std::min_element(mc_adc_values.cbegin(), mc_adc_values.cend());
402 *std::max_element(mc_adc_values.cbegin(), mc_adc_values.cend());
407 int mc_n_hits_in_slice4d = 0;
408 int mc_n_hits_in_slice4d_ge_100_ADC = 0;
411 for (
auto const& nova_slice : *nova_slices)
413 if (!nova_slice.IsNoise())
415 for (
auto const& nova_slice_hit : nova_slice.AllCells())
417 if (nova_slice_hit->IsMC() &&
420 ++mc_n_hits_in_slice4d;
422 if (nova_slice_hit->ADC() >= 100)
423 ++mc_n_hits_in_slice4d_ge_100_ADC;
431 int mc_n_hits_on_cosmic_track = 0;
432 int mc_n_hits_on_cosmic_track_ge_100_ADC = 0;
435 for (
auto const& cosmic_track : *cosmic_tracks)
437 for (
auto const& cosmic_track_hit : cosmic_track.AllCells())
439 if (cosmic_track_hit->IsMC() &&
442 ++mc_n_hits_on_cosmic_track;
444 if (cosmic_track_hit->ADC() >= 100)
445 ++mc_n_hits_on_cosmic_track_ge_100_ADC;
463 tset(
"event_length", event_length);
467 tset(
"ddt_trigger_decisions", ddt_decisions->size());
470 tset(
"n_hits_ge_100_ADC", n_hits_ge_100_ADC);
471 tset(
"n_hits_ge_150_ADC", n_hits_ge_150_ADC);
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);
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);
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);
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);
493 if (!offline_hits.failedToGet())
494 tset(
"n_hits", offline_hits->size());
498 tset(
"n_clusters", mono_clusters->size());
500 if (mono_clusters->size() == 1)
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());
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());
522 tset(
"n_slices_higheremoved",
n_slices(*higheremoved_slices));
526 find_hough_results(higheremoved_slices, e,
"hough");
529 std::map<geo::View_t, unsigned> n_peaks_max =
533 for (
size_t n_slice = 0; n_slice != higheremoved_slices->size(); ++n_slice)
535 if (!higheremoved_slices->at(n_slice).IsNoise())
537 std::vector<art::Ptr<rb::HoughResult> > hough_results;
538 find_hough_results.
get(n_slice, hough_results);
540 for (
auto const& hough_result : hough_results)
542 auto view = hough_result->fView;
543 unsigned n_peak = hough_result->fPeak.size();
554 tset(
"reco_n_tracks", reco_tracks->size());
558 for (
unsigned n_track = 1; n_track <= n_tracks; ++n_track)
560 auto track = reco_tracks->at(n_track - 1);
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());
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());
586 tset(name +
"sum_adc",
track.cluster().TotalADC());
588 tset(name +
"tstart",
track.cluster().MinTNS());
589 tset(name +
"tend",
track.cluster().MaxTNS());
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());
604 tset(name +
"dt",
track.cluster().ExtentTNS());
606 tset(name +
"max_time_gap_xz",
track.max_time_gap_xz());
607 tset(name +
"max_time_gap_yz",
track.max_time_gap_yz());
609 tset(name +
"most_missing_planes_xz",
611 tset(name +
"most_missing_planes_yz",
619 tset(
"mc_monopole_hit_detector",
false);
641 if (gen_info->size() == 1)
644 unsigned const trajectory = 0;
646 auto m = particle.
Mass();
649 tset(
"gen_dxdz",
p.Px() /
p.Pz());
650 tset(
"gen_dydz",
p.Py() /
p.Pz());
651 tset(
"gen_beta",
p.P() /
m);
654 << __FILE__ <<
":" << __LINE__ <<
"\n" 655 <<
"The generator object contains " << gen_info->size()
656 <<
" elements, there should only be 1." <<
std::endl;
661 std::cout <<
"MF: Number of MC Particles in Navigator = " 665 tset(
"mc_monopole_hit_detector",
true);
668 unsigned const trajectory = 0;
670 auto m = particle->
Mass();
676 tset(
"mc_dxdz",
p.Px() /
p.Pz());
677 tset(
"mc_dydz",
p.Py() /
p.Pz());
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);
685 std::map<std::string, std::vector<sim::FLSHit> > fls_hits;
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"])
696 double path_length =
hit.GetTotalPathLength();
697 if (path_length > 5.0)
699 Edep_sum +=
hit.GetEdep();
700 Edep_Birks_sum +=
hit.GetEdepBirks();
701 path_length_sum += path_length;
703 h_.at(
"dEdx")->Fill(
hit.GetEdep() / path_length * 1e3);
707 tset(
"mc_path_length_sum", path_length_sum);
708 if (path_length_sum > 0)
710 tset(
"mc_dEdx", Edep_sum / path_length_sum);
711 tset(
"mc_dEdx_Birks", Edep_Birks_sum / path_length_sum);
712 tset(
"mc_n_cerenkov", n_cerenkov_sum / path_length_sum);
715 for (
auto const&
hits : fls_hits)
719 tset(
"mc_n_hits_" + suffix,
hits.second.size());
724 std::sort(fls_hits.at(
"all").begin(), fls_hits.at(
"all").end(),
726 {
return lhs.
GetEntryT() < rhs.GetEntryT(); });
728 auto entry_hit = fls_hits.at(
"all").front();
729 auto exit_hit = fls_hits.at(
"all").back();
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());
740 std::map<std::string, std::vector<double> > edeps;
741 for (
auto const&
hit : offline_hitptrs)
744 for (
auto const& fls_hit : fls_hits)
747 static_cast<double>(
hit->ADC()) / fls_hit.GetTotalPathLength();
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);
765 double average_edep = sum /
static_cast<double>(
position.second.size());
766 tset(
"mc_" +
position.first +
"_average_edep", average_edep);
770 tset(
"is_mc",
false);
779 (std::vector<sim::FLSHit>
const&
hits)
const 784 auto it = std::minmax_element
785 (hits.cbegin(), hits.cend(),
786 [](
auto const lhs,
auto const rhs)
787 {
return lhs.GetCellID() < rhs.GetCellID(); });
790 return std::abs(
it.first->GetCellID() -
it.second->GetCellID()) + 1;
796 (std::vector<sim::FLSHit>
const& hits)
const 801 auto it = std::minmax_element
802 (hits.cbegin(), hits.cend(),
803 [](
auto const lhs,
auto const rhs)
804 {
return lhs.GetPlaneID() < rhs.GetPlaneID(); });
807 return std::abs(
it.first->GetPlaneID() -
it.second->GetPlaneID()) + 1;
841 unsigned result = std::count_if
842 (slices.begin(), slices.end(),
843 [](
rb::Cluster const& slice){
return !slice.IsNoise(); });
852 double Px(particle->
Px()), Py(particle->
Py());
865 for (
auto const &
hit : hitmap.at(
"all"))
868 hitmap.at(
"x").push_back(
hit);
870 hitmap.at(
"y").push_back(
hit);
878 double Px(particle->
Px()), Py(particle->
Py()), Pz(particle->
Pz());
879 double Pt =
std::sqrt(Px * Px + Py * Py);
T max(const caf::Proxy< T > &a, T b)
SubRunNumber_t subRun() const
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
int GetPlaneID() const
Plane ID.
const sim::ParticleNavigator & ParticleNavigator() const
Get a reference to the ParticleNavigator.
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
art::ServiceHandle< art::TFileService > tfs_
int GetCellID() const
Cell ID.
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...
TVector3 get_position(sim::FLSHit const &hit) const
double Px(const int i=0) const
unsigned n_slices(std::vector< rb::Cluster > const &slices) const
std::string const & fileName() const
A collection of associated CellHits.
A single unit of energy deposition in the liquid scintillator.
const double SPEED_OF_LIGHT
const PlaneGeo * Plane(unsigned int i) const
DEFINE_ART_MODULE(TestTMapFile)
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.
bool IsNoise(const art::Ptr< rb::CellHit > &hit) const
Is this hit not associated with any particles?
int extent_plane(std::vector< sim::FLSHit > const &hits) const
double branch_default_value_
const sim::Particle * Primary(const int) const
const XML_Char int const XML_Char * value
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
Monopole & operator=(Monopole const &)=delete
double phi(sim::Particle const *particle) const
void analyze(art::Event const &e) override
EventNumber_t event() const
EDAnalyzer(Table< Config > const &config)
unsigned n_mc_hits(mono::Track3D const &track) const
void respondToOpenInputFile(art::FileBlock const &fb) override
cet::coded_exception< errors::ErrorCodes, ExceptionDetail::translate > Exception
std::map< std::string, double > t_
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
std::map< std::string, TH1 * > h_
double pythag(double x, double y)
2D Euclidean distance
const TLorentzVector & Momentum(const int i=0) const
bool HaveTruthInfo() const
Is this a file with truth info in? (Is BackTracker going to be any use to you?)
double Pz(const int i=0) const
Monopole(fhicl::ParameterSet const &p)
size_type get(size_type i, reference item, data_reference data) const
geo::View_t view(unsigned plane) const
std::string to_string(ModuleType mt)
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.
art::ServiceHandle< geo::Geometry > geo_
art::ServiceHandle< cheat::BackTracker > bt_
Encapsulate the geometry of one entire detector (near, far, ndos)
int extent_cell(std::vector< sim::FLSHit > const &hits) const