All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Public Member Functions | Private Member Functions | Private Attributes | List of all members
numu::MuPVertexStudy Class Reference
Inheritance diagram for numu::MuPVertexStudy:

Public Member Functions

 MuPVertexStudy (fhicl::ParameterSet const &p)
 
 MuPVertexStudy (MuPVertexStudy const &)=delete
 
 MuPVertexStudy (MuPVertexStudy &&)=delete
 
MuPVertexStudy & operator= (MuPVertexStudy const &)=delete
 
MuPVertexStudy & operator= (MuPVertexStudy &&)=delete
 
void analyze (art::Event const &e) override
 

Private Member Functions

void InitTrue ()
 
void InitReco (unsigned i, const std::string &prefix)
 
void InitVertex ()
 
void FillReco (unsigned i, const simb::MCParticle &muon, const simb::MCParticle &proton, float muon_completion, float proton_completion, float muon_purity, float proton_purity)
 
void FillVertex (const TVector3 vertex, const recob::Track *muon, const recob::Track *proton)
 
void FillTrue (const simb::MCParticle &muon, const simb::MCParticle &proton)
 
float ParticleLength (const simb::MCParticle &particle)
 

Private Attributes

TrueHistos fTrueHistos
 
std::array< TrackHistos, 3 > fTrackHistos
 
VertexHistos fVertexHistos
 
std::vector< std::vector
< geo::BoxBoundedGeo > > 		fTPCVolumes
 
std::vector< geo::BoxBoundedGeo > fActiveVolumes
 

Detailed Description

Definition at line 86 of file MuPVertexStudy_module.cc.

Constructor & Destructor Documentation

numu::MuPVertexStudy::MuPVertexStudy ( fhicl::ParameterSet const &  p)
explicit

Definition at line 148 of file MuPVertexStudy_module.cc.

149  : EDAnalyzer{p} // ,
150  // More initializers here.
151 {
152 
153  InitTrue();
154  InitVertex();
155 
156  std::vector<std::string> prefixes {"both", "mu", "p"};
157  for (unsigned i = 0; i < 3; i++) {
158  InitReco(i, prefixes[i]);
159  }
160 
161  const geo::GeometryCore *geometry = lar::providerFrom<geo::Geometry>();
162 
163  // first the TPC volumes
164  for (auto const &cryo: geometry->IterateCryostats()) {
165  geo::GeometryCore::TPC_iterator iTPC = geometry->begin_TPC(cryo.ID()),
166  tend = geometry->end_TPC(cryo.ID());
167  std::vector<geo::BoxBoundedGeo> this_tpc_volumes;
168  while (iTPC != tend) {
169  geo::TPCGeo const& TPC = *iTPC;
170  this_tpc_volumes.push_back(TPC.ActiveBoundingBox());
171  iTPC++;
172  }
173  fTPCVolumes.push_back(std::move(this_tpc_volumes));
174  }
175 
176  // then combine them into active volumes
177  for (const std::vector<geo::BoxBoundedGeo> &tpcs: fTPCVolumes) {
178  double XMin = std::min_element(tpcs.begin(), tpcs.end(), [](auto &lhs, auto &rhs) { return lhs.MinX() < rhs.MinX(); })->MinX();
179  double YMin = std::min_element(tpcs.begin(), tpcs.end(), [](auto &lhs, auto &rhs) { return lhs.MinY() < rhs.MinY(); })->MinY();
180  double ZMin = std::min_element(tpcs.begin(), tpcs.end(), [](auto &lhs, auto &rhs) { return lhs.MinZ() < rhs.MinZ(); })->MinZ();
181 
182  double XMax = std::max_element(tpcs.begin(), tpcs.end(), [](auto &lhs, auto &rhs) { return lhs.MaxX() < rhs.MaxX(); })->MaxX();
183  double YMax = std::max_element(tpcs.begin(), tpcs.end(), [](auto &lhs, auto &rhs) { return lhs.MaxY() < rhs.MaxY(); })->MaxY();
184  double ZMax = std::max_element(tpcs.begin(), tpcs.end(), [](auto &lhs, auto &rhs) { return lhs.MaxZ() < rhs.MaxZ(); })->MaxZ();
185 
186  fActiveVolumes.emplace_back(XMin, XMax, YMin, YMax, ZMin, ZMax);
187  }
188 
189 }
numu::MuPVertexStudy::InitReco
void InitReco(unsigned i, const std::string &prefix)
Definition: MuPVertexStudy_module.cc:129
geometry
const geo::GeometryCore * geometry
Definition: NuVertexChargeTree_module.cc:373
p
pdgs p
Definition: selectors.fcl:22
geo::TPCGeo
Geometry information for a single TPC.
Definition: TPCGeo.h:38
numu::MuPVertexStudy::fActiveVolumes
std::vector< geo::BoxBoundedGeo > fActiveVolumes
Definition: MuPVertexStudy_module.cc:107
geo::TPCGeo::ActiveBoundingBox
geo::BoxBoundedGeo const & ActiveBoundingBox() const
Returns the box of the active volume of this TPC.
Definition: TPCGeo.h:320
TPC
BEGIN_PROLOG TPC
Definition: sbnd_flashalgo.fcl:14
geo::GeometryCore::begin_TPC
TPC_iterator begin_TPC() const
Returns an iterator pointing to the first TPC in the detector.
Definition: GeometryCore.h:2570
numu::MuPVertexStudy::fTPCVolumes
std::vector< std::vector< geo::BoxBoundedGeo > > fTPCVolumes
Definition: MuPVertexStudy_module.cc:106
geo::GeometryCore
Description of geometry of one entire detector.
Definition: GeometryCore.h:1476
geo::GeometryCore::IterateCryostats
IteratorBox< cryostat_iterator,&GeometryCore::begin_cryostat,&GeometryCore::end_cryostat > IterateCryostats() const
Enables ranged-for loops on all cryostats of the detector.
Definition: GeometryCore.h:2206
geo::details::geometry_element_iterator
Forward iterator browsing all geometry elements in the detector.
Definition: GeometryCore.h:727
geo::GeometryCore::end_TPC
TPC_iterator end_TPC() const
Returns an iterator pointing after the last TPC in the detector.
Definition: GeometryCore.h:2574
numu::MuPVertexStudy::MuPVertexStudy ( MuPVertexStudy const &  )
delete
numu::MuPVertexStudy::MuPVertexStudy ( MuPVertexStudy &&  )
delete

Member Function Documentation

void numu::MuPVertexStudy::analyze ( art::Event const &  e)
override

Definition at line 235 of file MuPVertexStudy_module.cc.

236 {
237  auto const clock_data = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(ev);
238  art::ServiceHandle<cheat::BackTrackerService> backtracker;
239 
240  const std::vector<simb::MCParticle> &mcparticle_list = *ev.getValidHandle<std::vector<simb::MCParticle>>("largeant");
241 
242  const simb::MCParticle *muon = NULL;
243  const simb::MCParticle *proton = NULL;
244 
245  // get the muon and proton
246  for (unsigned i = 0; i < mcparticle_list.size(); i++) {
247  if (mcparticle_list[i].PdgCode() == 13 && mcparticle_list[i].Process() == "primary") {
248  assert(muon == NULL);
249  muon = &mcparticle_list[i];
250  }
251  if (mcparticle_list[i].PdgCode() == 2212 && mcparticle_list[i].Process() == "primary") {
252  assert(proton == NULL);
253  proton = &mcparticle_list[i];
254  }
255  }
256 
257  assert(muon != NULL);
258  assert(proton != NULL);
259 
260  FillTrue(*muon, *proton);
261 
262  float muonE = 0.;
263  std::vector<const sim::IDE*> muon_ides(backtracker->TrackIdToSimIDEs_Ps(muon->TrackId()));
264  for (const sim::IDE *ide: muon_ides) muonE += ide->energy;
265 
266  float protonE = 0.;
267  std::vector<const sim::IDE*> proton_ides(backtracker->TrackIdToSimIDEs_Ps(proton->TrackId()));
268  for (const sim::IDE *ide: proton_ides) protonE += ide->energy;
269 
270  const auto &particle_handle = ev.getValidHandle<std::vector<recob::PFParticle>>("pandora");
271  const std::vector<recob::PFParticle> &particle_list = *particle_handle;
272  art::FindManyP<recob::Cluster> particles_to_clusters(particle_handle, ev, "pandora");
273  art::FindManyP<recob::Track> particles_to_tracks(particle_handle, ev, "pandoraTrack");
274 
275  const recob::Track *muon_track = NULL;
276  const recob::Track *proton_track = NULL;
277 
278  for (unsigned i_part = 0; i_part < particle_list.size(); i_part++) {
279  const std::vector<art::Ptr<recob::Cluster>> clusters = particles_to_clusters.at(i_part);
280  art::FindManyP<recob::Hit> clusters_to_hits(clusters, ev, "pandora");
281 
282  std::vector<art::Ptr<recob::Hit>> hits;
283  for (unsigned i_clus = 0; i_clus < clusters.size(); i_clus++) {
284  const std::vector<art::Ptr<recob::Hit>> &this_hits = clusters_to_hits.at(i_clus);
285  hits.insert(hits.end(), this_hits.begin(), this_hits.end());
286  }
287 
288  std::vector<std::pair<int, float>> matches = CAFRecoUtils::AllTrueParticleIDEnergyMatches(clock_data, hits);
289  float totalE = CAFRecoUtils::TotalHitEnergy(clock_data, hits);
290 
291  float muon_completion = Completion(*muon, muonE, matches, mcparticle_list);
292  float proton_completion = Completion(*proton, protonE, matches, mcparticle_list);
293 
294  float muon_purity = Purity(*muon, totalE, matches, mcparticle_list);
295  float proton_purity = Purity(*proton, totalE, matches, mcparticle_list);
296 
297  const recob::Track *track = particles_to_tracks.at(i_part).size() ? particles_to_tracks.at(i_part).at(0).get() : NULL;
298 
299  if (muon_completion > 0.5 && proton_completion > 0.5) {
300  FillReco(0, *muon, *proton, muon_completion, proton_completion, muon_purity, proton_purity);
301  muon_track = track;
302  proton_track = track;
303  }
304  else if (muon_completion > 0.5) {
305  FillReco(1, *muon, *proton, muon_completion, proton_completion, muon_purity, proton_purity);
306  muon_track = track;
307  }
308  else if (proton_completion > 0.5) {
309  FillReco(2, *muon, *proton, muon_completion, proton_completion, muon_purity, proton_purity);
310  proton_track = track;
311  }
312  }
313  if (muon_track && proton_track && muon_track != proton_track) {
314  FillVertex(muon->Position().Vect(), muon_track, proton_track);
315  }
316 }
numu::MuPVertexStudy::FillTrue
void FillTrue(const simb::MCParticle &muon, const simb::MCParticle &proton)
Definition: MuPVertexStudy_module.cc:208
numu::MuPVertexStudy::FillReco
void FillReco(unsigned i, const simb::MCParticle &muon, const simb::MCParticle &proton, float muon_completion, float proton_completion, float muon_purity, float proton_purity)
Definition: MuPVertexStudy_module.cc:215
track
process_name use argoneut_mc_hitfinder track
Definition: standard_reco.fcl:36
CAFRecoUtils::AllTrueParticleIDEnergyMatches
std::vector< std::pair< int, float > > AllTrueParticleIDEnergyMatches(const detinfo::DetectorClocksData &clockData, const std::vector< art::Ptr< recob::Hit > > &hits, bool rollup_unsaved_ids=1)
Definition: sbncode/sbncode/CAFMaker/RecoUtils/RecoUtils.cc:3
sim::IDE
Ionization at a point of the TPC sensitive volume.
Definition: SimChannel.h:86
muon
process_name can override from command line with o or output muon
Definition: runPID.fcl:28
Completion
float Completion(const simb::MCParticle &particle, float particleE, const std::vector< std::pair< int, float >> &matches, const std::vector< simb::MCParticle > &particles)
Definition: MuPVertexStudy_module.cc:318
Purity
float Purity(const simb::MCParticle &particle, float totalE, const std::vector< std::pair< int, float >> &matches, const std::vector< simb::MCParticle > &particles)
Definition: MuPVertexStudy_module.cc:327
numu::MuPVertexStudy::FillVertex
void FillVertex(const TVector3 vertex, const recob::Track *muon, const recob::Track *proton)
Definition: MuPVertexStudy_module.cc:229
CAFRecoUtils::TotalHitEnergy
float TotalHitEnergy(const detinfo::DetectorClocksData &clockData, const std::vector< art::Ptr< recob::Hit > > &hits)
Definition: sbncode/sbncode/CAFMaker/RecoUtils/RecoUtils.cc:47
recob::Track
Track from a non-cascading particle.A recob::Track consists of a recob::TrackTrajectory, plus additional members relevant for a &quot;fitted&quot; track:
Definition: bj/lardataobj/RecoBase/Track.h:49
void numu::MuPVertexStudy::FillReco ( unsigned  i,
const simb::MCParticle &  muon,
const simb::MCParticle &  proton,
float  muon_completion,
float  proton_completion,
float  muon_purity,
float  proton_purity 
)
private

Definition at line 215 of file MuPVertexStudy_module.cc.

215  {
216 
217  fTrackHistos[i].muon_completion->Fill(muon_completion);
218  fTrackHistos[i].proton_completion->Fill(proton_completion);
219  fTrackHistos[i].muon_purity->Fill(muon_purity);
220  fTrackHistos[i].proton_purity->Fill(proton_purity);
221 
222  fTrackHistos[i].muon_length->Fill(ParticleLength(muon));
223  fTrackHistos[i].proton_length->Fill(ParticleLength(proton));
224 
225  fTrackHistos[i].cos_open_angle->Fill(muon.Momentum().Vect().Dot(proton.Momentum().Vect()) / (muon.Momentum().Vect().Mag() * proton.Momentum().Vect().Mag()));
226 
227 }
numu::MuPVertexStudy::ParticleLength
float ParticleLength(const simb::MCParticle &particle)
Definition: MuPVertexStudy_module.cc:191
muon
process_name can override from command line with o or output muon
Definition: runPID.fcl:28
numu::MuPVertexStudy::fTrackHistos
std::array< TrackHistos, 3 > fTrackHistos
Definition: MuPVertexStudy_module.cc:104
void numu::MuPVertexStudy::FillTrue ( const simb::MCParticle &  muon,
const simb::MCParticle &  proton 
)
private

Definition at line 208 of file MuPVertexStudy_module.cc.

208  {
209  fTrueHistos.muon_length->Fill(ParticleLength(muon));
210  fTrueHistos.proton_length->Fill(ParticleLength(proton));
211 
212  fTrueHistos.cos_open_angle->Fill(muon.Momentum().Vect().Dot(proton.Momentum().Vect()) / (muon.Momentum().Vect().Mag() * proton.Momentum().Vect().Mag()));
213 }
numu::MuPVertexStudy::ParticleLength
float ParticleLength(const simb::MCParticle &particle)
Definition: MuPVertexStudy_module.cc:191
muon
process_name can override from command line with o or output muon
Definition: runPID.fcl:28
void numu::MuPVertexStudy::FillVertex ( const TVector3  vertex,
const recob::Track *  muon,
const recob::Track *  proton 
)
private

Definition at line 229 of file MuPVertexStudy_module.cc.

229  {
230  TVector3 muon_start(muon->Start().X(), muon->Start().Y(), muon->Start().Z());
231  TVector3 muon_end(muon->End().X(), muon->End().Y(), muon->End().Z());
232  fVertexHistos.vertex_diff->Fill(std::min((vertex - muon_start).Mag(), (vertex - muon_end).Mag()));
233 }
vertex
process_name vertex
Definition: cheaterreco.fcl:51
recob::Track::Start
Point_t const & Start() const
Access to track position at different points.
Definition: bj/lardataobj/RecoBase/Track.h:123
recob::Track::End
Point_t const & End() const
Definition: bj/lardataobj/RecoBase/Track.h:125
void numu::MuPVertexStudy::InitReco ( unsigned  i,
const std::string &  prefix 
)
private

Definition at line 129 of file MuPVertexStudy_module.cc.

129  {
130  art::ServiceHandle<art::TFileService> tfs;
131 
132  fTrackHistos[i].cos_open_angle = tfs->make<TH1D>((prefix + "reco_CosOpenAngle").c_str(), "", 64, -3.2, 3.2);
133  fTrackHistos[i].muon_length = tfs->make<TH1D>((prefix + "reco_mu_length").c_str(), "", 100, 0., 500.);
134  fTrackHistos[i].proton_length = tfs->make<TH1D>((prefix + "reco_p_length").c_str(), "", 100, 0., 200.);
135 
136  fTrackHistos[i].muon_completion = tfs->make<TH1D>((prefix + "reco_mu_completion").c_str(), "", 100, 0., 1.);
137  fTrackHistos[i].proton_completion = tfs->make<TH1D>((prefix + "reco_p_completion").c_str(), "", 100, 0., 1.);
138  fTrackHistos[i].muon_purity = tfs->make<TH1D>((prefix + "reco_mu_purity").c_str(), "", 100, 0., 1.);
139  fTrackHistos[i].proton_purity = tfs->make<TH1D>((prefix + "reco_p_purity").c_str(), "", 100, 0., 1.);
140 }
numu::MuPVertexStudy::fTrackHistos
std::array< TrackHistos, 3 > fTrackHistos
Definition: MuPVertexStudy_module.cc:104
tfs
art::ServiceHandle< art::TFileService > tfs
Definition: NuVertexChargeTree_module.cc:254
void numu::MuPVertexStudy::InitTrue ( )
private

Definition at line 121 of file MuPVertexStudy_module.cc.

121  {
122  art::ServiceHandle<art::TFileService> tfs;
123 
124  fTrueHistos.cos_open_angle = tfs->make<TH1D>("CosOpenAngle", "", 64, -3.2, 3.2);
125  fTrueHistos.muon_length = tfs->make<TH1D>("mu_length", "", 100, 0., 500.);
126  fTrueHistos.proton_length = tfs->make<TH1D>("p_length", "", 100, 0., 200.);
127 }
tfs
art::ServiceHandle< art::TFileService > tfs
Definition: NuVertexChargeTree_module.cc:254
void numu::MuPVertexStudy::InitVertex ( )
private

Definition at line 142 of file MuPVertexStudy_module.cc.

142  {
143  art::ServiceHandle<art::TFileService> tfs;
144 
145  fVertexHistos.vertex_diff = tfs->make<TH1D>("vertex_diff", "", 200, 0, 100.);
146 }
tfs
art::ServiceHandle< art::TFileService > tfs
Definition: NuVertexChargeTree_module.cc:254
MuPVertexStudy& numu::MuPVertexStudy::operator= ( MuPVertexStudy const &  )
delete
MuPVertexStudy& numu::MuPVertexStudy::operator= ( MuPVertexStudy &&  )
delete
float numu::MuPVertexStudy::ParticleLength ( const simb::MCParticle &  particle)
private

Definition at line 191 of file MuPVertexStudy_module.cc.

191  {
192  std::vector<geoalgo::AABox> aa_volumes;
193  for (const geo::BoxBoundedGeo &AV: fActiveVolumes) {
194  if (AV.ContainsPosition(particle.Position().Vect())) {
195  aa_volumes.emplace_back(AV.MinX(), AV.MinY(), AV.MinZ(), AV.MaxX(), AV.MaxY(), AV.MaxZ());
196  break;
197  }
198  }
199 
200  float length = 0.;
201  for (unsigned i = 1; i < particle.NumberTrajectoryPoints(); i++) {
202  length += ContainedLength(particle.Position(i-1).Vect(), particle.Position(i).Vect(), aa_volumes);
203  }
204 
205  return length;
206 }
util::ContainedLength
double ContainedLength(const TVector3 &v0, const TVector3 &v1, const std::vector< geo::BoxBoundedGeo > &boxes)
Definition: Interaction.cxx:28
numu::MuPVertexStudy::fActiveVolumes
std::vector< geo::BoxBoundedGeo > fActiveVolumes
Definition: MuPVertexStudy_module.cc:107
geo::BoxBoundedGeo
A base class aware of world box coordinatesAn object describing a simple shape can inherit from this ...
Definition: BoxBoundedGeo.h:33

Member Data Documentation

std::vector<geo::BoxBoundedGeo> numu::MuPVertexStudy::fActiveVolumes
private

Definition at line 107 of file MuPVertexStudy_module.cc.

std::vector<std::vector<geo::BoxBoundedGeo> > numu::MuPVertexStudy::fTPCVolumes
private

Definition at line 106 of file MuPVertexStudy_module.cc.

std::array<TrackHistos, 3> numu::MuPVertexStudy::fTrackHistos
private

Definition at line 104 of file MuPVertexStudy_module.cc.

TrueHistos numu::MuPVertexStudy::fTrueHistos
private

Definition at line 103 of file MuPVertexStudy_module.cc.

VertexHistos numu::MuPVertexStudy::fVertexHistos
private

Definition at line 105 of file MuPVertexStudy_module.cc.

The documentation for this class was generated from the following file:
Generated by   doxygen 1.8.5