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
DUNE::NeutrinoTrackingEff Class Reference
Inheritance diagram for DUNE::NeutrinoTrackingEff:

Public Member Functions

 NeutrinoTrackingEff (fhicl::ParameterSet const &pset)
 

Private Member Functions

void beginJob ()
 
void endJob ()
 
void beginRun (const art::Run &run)
 
void analyze (const art::Event &evt)
 
void processEff (const art::Event &evt)
 
void truthMatcher (detinfo::DetectorClocksData const &clockData, std::vector< art::Ptr< recob::Hit >> all_hits, std::vector< art::Ptr< recob::Hit >> track_hits, const simb::MCParticle *&MCparticle, double &Efrac, double &Ecomplet)
 
double truthLength (const detinfo::DetectorClocksData &clockData, detinfo::DetectorPropertiesData const &detProp, const simb::MCParticle *MCparticle)
 
bool insideFV (double vertex[4])
 
void doEfficiencies ()
 

Private Attributes

std::string fMCTruthModuleLabel
 
std::string fTrackModuleLabel
 
int fNeutrinoPDGcode
 
int fLeptonPDGcode
 
double fMaxNeutrinoE
 
double fMaxLeptonP
 
bool fisNeutrinoInt
 
int MC_isCC
 
int MC_incoming_PDG
 
double MC_incoming_P [4]
 
double MC_vertex [4]
 
double MC_lepton_startMomentum [4]
 
int MC_leading_protonID
 
int MC_leading_PionPlusID
 
int MC_leading_PionMinusID
 
int MC_leptonID
 
int MC_kaonID
 
int MC_michelID
 
double MC_leptonP
 
double MC_leading_PionPlusP
 
double MC_leading_ProtonP
 
double MC_leading_PionMinusP
 
double MC_kaonP
 
double MC_michelP
 
TH1D * h_Ev_den
 
TH1D * h_Ev_num
 
TH1D * h_Pmu_den
 
TH1D * h_Pmu_num
 
TH1D * h_theta_den
 
TH1D * h_theta_num
 
TH1D * h_Pproton_den
 
TH1D * h_Pproton_num
 
TH1D * h_Ppion_plus_den
 
TH1D * h_Ppion_plus_num
 
TH1D * h_Ppion_minus_den
 
TH1D * h_Ppion_minus_num
 
TH1D * h_Efrac_lepton
 
TH1D * h_Ecomplet_lepton
 
TH1D * h_Efrac_proton
 
TH1D * h_Ecomplet_proton
 
TH1D * h_Efrac_pion_plus
 
TH1D * h_Ecomplet_pion_plus
 
TH1D * h_Efrac_pion_minus
 
TH1D * h_Ecomplet_pion_minus
 
TH1D * h_trackRes_lepton
 
TH1D * h_trackRes_proton
 
TH1D * h_trackRes_pion_plus
 
TH1D * h_trackRes_pion_minus
 
TH1D * h_muon_length
 
TH1D * h_proton_length
 
TH1D * h_pionp_length
 
TH1D * h_pionm_length
 
TH1D * h_muonwtrk_length
 
TH1D * h_protonwtrk_length
 
TH1D * h_pionpwtrk_length
 
TH1D * h_pionmwtrk_length
 
TEfficiency * h_Eff_Ev = 0
 
TEfficiency * h_Eff_Pmu = 0
 
TEfficiency * h_Eff_theta = 0
 
TEfficiency * h_Eff_Pproton = 0
 
TEfficiency * h_Eff_Ppion_plus = 0
 
TEfficiency * h_Eff_Ppion_minus = 0
 
TEfficiency * h_Eff_Lmuon = 0
 
TEfficiency * h_Eff_Lproton = 0
 
TEfficiency * h_Eff_Lpion_plus = 0
 
TEfficiency * h_Eff_Lpion_minus = 0
 
TH1D * h_Pkaon_den
 
TH1D * h_Pkaon_num
 
TH1D * h_Pmichel_e_den
 
TH1D * h_Pmichel_e_num
 
TH1D * h_Efrac_kaon
 
TH1D * h_Ecomplet_kaon
 
TH1D * h_trackRes_kaon
 
TH1D * h_Efrac_michel
 
TH1D * h_Ecomplet_michel
 
TH1D * h_trackRes_michel
 
TH1D * h_kaon_length
 
TH1D * h_michel_length
 
TH1D * h_kaonwtrk_length
 
TH1D * h_michelwtrk_length
 
TEfficiency * h_Eff_Pkaon = 0
 
TEfficiency * h_Eff_Pmichel = 0
 
TEfficiency * h_Eff_Lkaon = 0
 
TEfficiency * h_Eff_Lmichel = 0
 
float fFidVolCutX
 
float fFidVolCutY
 
float fFidVolCutZ
 
float fFidVolXmin
 
float fFidVolXmax
 
float fFidVolYmin
 
float fFidVolYmax
 
float fFidVolZmin
 
float fFidVolZmax
 
double fDriftVelocity
 
art::ServiceHandle
< geo::Geometry const > 
geom
 

Detailed Description

Definition at line 51 of file NeutrinoTrackingEff_module.cc.

Constructor & Destructor Documentation

DUNE::NeutrinoTrackingEff::NeutrinoTrackingEff ( fhicl::ParameterSet const &  pset)
explicit

Definition at line 187 of file NeutrinoTrackingEff_module.cc.

187  : EDAnalyzer(p)
188  {
189  fMCTruthModuleLabel = p.get<std::string>("MCTruthModuleLabel");
190  fTrackModuleLabel = p.get<std::string>("TrackModuleLabel");
191  fisNeutrinoInt = p.get<bool>("isNeutrinoInt");
192  fLeptonPDGcode = p.get<int>("LeptonPDGcode");
193  fNeutrinoPDGcode = p.get<int>("NeutrinoPDGcode");
194  fMaxNeutrinoE = p.get<double>("MaxNeutrinoE");
195  fMaxLeptonP = p.get<double>("MaxLeptonP");
196  fFidVolCutX = p.get<float>("FidVolCutX");
197  fFidVolCutY = p.get<float>("FidVolCutY");
198  fFidVolCutZ = p.get<float>("FidVolCutZ");
199 
200  auto const detProp =
201  art::ServiceHandle<detinfo::DetectorPropertiesService const>()->DataForJob();
202  fDriftVelocity = detProp.DriftVelocity() * 1e-3; // cm/ns
203  }
pdgs p
Definition: selectors.fcl:22
do i e
auto const detProp

Member Function Documentation

void DUNE::NeutrinoTrackingEff::analyze ( const art::Event &  evt)
private

Definition at line 428 of file NeutrinoTrackingEff_module.cc.

429  {
430  if (event.isRealData()) return;
431 
432  processEff(event);
433  }
void processEff(const art::Event &evt)
void DUNE::NeutrinoTrackingEff::beginJob ( )
private

Definition at line 206 of file NeutrinoTrackingEff_module.cc.

207  {
208  std::cout << "job begin..." << std::endl;
209  auto const* geo = lar::providerFrom<geo::Geometry>();
210  // Define histogram boundaries (cm).
211  // For now only draw cryostat=0.
212  double minx = 1e9;
213  double maxx = -1e9;
214  double miny = 1e9;
215  double maxy = -1e9;
216  double minz = 1e9;
217  double maxz = -1e9;
218  for (size_t i = 0; i < geo->NTPC(); ++i) {
219  double local[3] = {0., 0., 0.};
220  double world[3] = {0., 0., 0.};
221  const geo::TPCGeo& tpc = geo->TPC(i);
222  tpc.LocalToWorld(local, world);
223  if (minx > world[0] - geo->DetHalfWidth(i)) minx = world[0] - geo->DetHalfWidth(i);
224  if (maxx < world[0] + geo->DetHalfWidth(i)) maxx = world[0] + geo->DetHalfWidth(i);
225  if (miny > world[1] - geo->DetHalfHeight(i)) miny = world[1] - geo->DetHalfHeight(i);
226  if (maxy < world[1] + geo->DetHalfHeight(i)) maxy = world[1] + geo->DetHalfHeight(i);
227  if (minz > world[2] - geo->DetLength(i) / 2.) minz = world[2] - geo->DetLength(i) / 2.;
228  if (maxz < world[2] + geo->DetLength(i) / 2.) maxz = world[2] + geo->DetLength(i) / 2.;
229  }
230 
231  fFidVolXmin = minx + fFidVolCutX;
232  fFidVolXmax = maxx - fFidVolCutX;
233  fFidVolYmin = miny + fFidVolCutY;
234  fFidVolYmax = maxy - fFidVolCutY;
235  fFidVolZmin = minz + fFidVolCutZ;
236  fFidVolZmax = maxz - fFidVolCutZ;
237 
238  std::cout << "Fiducial volume:"
239  << "\n"
240  << fFidVolXmin << "\t< x <\t" << fFidVolXmax << "\n"
241  << fFidVolYmin << "\t< y <\t" << fFidVolYmax << "\n"
242  << fFidVolZmin << "\t< z <\t" << fFidVolZmax << "\n";
243 
244  art::ServiceHandle<art::TFileService const> tfs;
245 
246  double E_bins[21] = {0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4, 4.5, 5.0,
247  5.5, 6.0, 7.0, 8.0, 10.0, 12.0, 14.0, 17.0, 20.0, 25.0};
248  double theta_bin[44] = {0., 1., 2., 3., 4., 5., 6., 7., 8., 9., 10.,
249  11., 12., 13., 14., 15., 16., 17., 18., 19., 20., 22.,
250  24., 26., 28., 30., 32., 34., 36., 38., 40., 42., 44.,
251  46., 48., 50., 55., 60., 65., 70., 75., 80., 85., 90.};
252  double Pbins[18] = {
253  0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.2, 1.4, 1.6, 1.8, 2.0, 2.5, 3.0};
254 
255  for (int i = 0; i < 21; ++i)
256  E_bins[i] *= fMaxNeutrinoE / 25.;
257  for (int i = 0; i < 18; ++i)
258  Pbins[i] *= fMaxLeptonP / 3.0;
259 
260  h_Ev_den = tfs->make<TH1D>(
261  "h_Ev_den", "Neutrino Energy; Neutrino Energy (GeV); Tracking Efficiency", 20, E_bins);
262  h_Ev_num = tfs->make<TH1D>(
263  "h_Ev_num", "Neutrino Energy; Neutrino Energy (GeV); Tracking Efficiency", 20, E_bins);
264  h_Pmu_den = tfs->make<TH1D>(
265  "h_Pmu_den", "Muon Momentum; Muon Momentum (GeV); Tracking Efficiency", 20, E_bins);
266  h_Pmu_num = tfs->make<TH1D>(
267  "h_Pmu_num", "Muon Momentum; Muon Momentum (GeV); Tracking Efficiency", 20, E_bins);
268  h_theta_den =
269  tfs->make<TH1D>("h_theta_den",
270  "Theta; Theta w.r.t beam direction (Degrees); Tracking Efficiency",
271  43,
272  theta_bin);
273  h_theta_num =
274  tfs->make<TH1D>("h_theta_num",
275  "Theta; Theta w.r.t beam direction (Degrees); Tracking Efficiency",
276  43,
277  theta_bin);
278  h_Pproton_den = tfs->make<TH1D>(
279  "h_Pproton_den", "Protons; Proton Momentum (GeV); Tracking Efficiency", 17, Pbins);
280  h_Pproton_num = tfs->make<TH1D>(
281  "h_Pproton_num", "Protons; Proton Momentum (GeV); Tracking Efficiency", 17, Pbins);
282  h_Ppion_plus_den = tfs->make<TH1D>(
283  "h_Ppion_plus_den", "Pions Plus; Pion Momentum (GeV); Tracking Efficiency", 17, Pbins);
284  h_Ppion_plus_num = tfs->make<TH1D>(
285  "h_Ppion_plus_num", "Pions Plus; Pion Momentum (GeV); Tracking Efficiency", 17, Pbins);
286  h_Ppion_minus_den = tfs->make<TH1D>(
287  "h_Ppion_minus_den", "Pions Minus; Pion Momentum (GeV); Tracking Efficiency", 17, Pbins);
288  h_Ppion_minus_num = tfs->make<TH1D>(
289  "h_Ppion_minus_num", "Pions Minus; Pion Momentum (GeV); Tracking Efficiency", 17, Pbins);
290  h_Ev_den->Sumw2();
291  h_Ev_num->Sumw2();
292  h_Pmu_den->Sumw2();
293  h_Pmu_num->Sumw2();
294  h_theta_den->Sumw2();
295  h_theta_num->Sumw2();
296  h_Pproton_den->Sumw2();
297  h_Pproton_num->Sumw2();
298  h_Ppion_plus_den->Sumw2();
299  h_Ppion_plus_num->Sumw2();
300  h_Ppion_minus_den->Sumw2();
301  h_Ppion_minus_num->Sumw2();
302 
303  h_Efrac_lepton = tfs->make<TH1D>("h_Efrac_lepton", "Efrac Lepton; Track Purity;", 60, 0, 1.2);
305  tfs->make<TH1D>("h_Ecomplet_lepton", "Ecomplet Lepton; Track Completeness;", 60, 0, 1.2);
306  h_Efrac_proton = tfs->make<TH1D>("h_Efrac_proton", "Efrac Proton; Track Purity;", 60, 0, 1.2);
308  tfs->make<TH1D>("h_Ecomplet_proton", "Ecomplet Proton; Track Completeness;", 60, 0, 1.2);
310  tfs->make<TH1D>("h_Efrac_pion_plus", "Efrac Pion +; Track Purity;", 60, 0, 1.2);
312  tfs->make<TH1D>("h_Ecomplet_pion_plus", "Ecomplet Pion +; Track Completeness;", 60, 0, 1.2);
314  tfs->make<TH1D>("h_Efrac_pion_minus", "Efrac Pion -; Track Purity;", 60, 0, 1.2);
316  tfs->make<TH1D>("h_Ecomplet_pion_minus", "Ecomplet Pion -; Track Completeness;", 60, 0, 1.2);
317  h_trackRes_lepton = tfs->make<TH1D>(
318  "h_trackRes_lepton", "Muon Residual; Truth length - Reco length (cm);", 200, -100, 100);
319  h_trackRes_proton = tfs->make<TH1D>(
320  "h_trackRes_proton", "Proton Residual; Truth length - Reco length (cm);", 200, -100, 100);
321  h_trackRes_pion_plus = tfs->make<TH1D>(
322  "h_trackRes_pion_plus", "Pion + Residual; Truth length - Reco length (cm);", 200, -100, 100);
323  h_trackRes_pion_minus = tfs->make<TH1D>(
324  "h_trackRes_pion_minus", "Pion - Residual; Truth length - Reco length (cm);", 200, -100, 100);
325  h_Efrac_lepton->Sumw2();
326  h_Ecomplet_lepton->Sumw2();
327  h_Efrac_proton->Sumw2();
328  h_Ecomplet_proton->Sumw2();
329  h_Efrac_pion_plus->Sumw2();
330  h_Ecomplet_pion_plus->Sumw2();
331  h_Efrac_pion_minus->Sumw2();
332  h_Ecomplet_pion_minus->Sumw2();
333  h_trackRes_lepton->Sumw2();
334  h_trackRes_proton->Sumw2();
335  h_trackRes_pion_plus->Sumw2();
336  h_trackRes_pion_minus->Sumw2();
337 
338  h_muon_length =
339  tfs->make<TH1D>("h_muon_length", "Muon Length; Muon Truth Length (cm)", 40, 0, 100);
341  tfs->make<TH1D>("h_proton_length", "Proton Length; Proton Truth Length (cm)", 40, 0, 100);
343  tfs->make<TH1D>("h_pionp_length", "Pion + Length; Pion^{+} Truth Length (cm)", 40, 0, 100);
345  tfs->make<TH1D>("h_pionm_length", "Pion - Length; Pion^{-} Truth Length (cm)", 40, 0, 100);
346 
348  tfs->make<TH1D>("h_muonwtrk_length", "Muon Length; Muon Truth Length (cm)", 40, 0, 100);
350  tfs->make<TH1D>("h_protonwtrk_length", "Proton Length; Proton Truth Length (cm)", 40, 0, 100);
351  h_pionpwtrk_length = tfs->make<TH1D>(
352  "h_pionpwtrk_length", "Pion + Length; Pion^{+} Truth Length (cm)", 40, 0, 100);
353  h_pionmwtrk_length = tfs->make<TH1D>(
354  "h_pionmwtrk_length", "Pion - Length; Pion^{-} Truth Length (cm)", 40, 0, 100);
355 
356  h_muon_length->Sumw2();
357  h_muonwtrk_length->Sumw2();
358  h_proton_length->Sumw2();
359  h_protonwtrk_length->Sumw2();
360  h_pionp_length->Sumw2();
361  h_pionpwtrk_length->Sumw2();
362  h_pionm_length->Sumw2();
363  h_pionmwtrk_length->Sumw2();
364 
365  h_Pkaon_den =
366  tfs->make<TH1D>("h_Pkaon_den", "Kaon; Kaon Momentum (GeV); Tracking Efficiency", 17, Pbins);
367  h_Pkaon_num =
368  tfs->make<TH1D>("h_Pkaon_num", "Kaon; Kaon Momentum (GeV); Tracking Efficiency", 17, Pbins);
370  tfs->make<TH1D>("h_Pmichel_e_den",
371  "Michel Electron; Michele e Momentum (GeV); Tracking Efficiency",
372  17,
373  Pbins);
375  tfs->make<TH1D>("h_Pmichel_e_num",
376  "Michel Electron; Michele e Momentum (GeV); Tracking Efficiency",
377  17,
378  Pbins);
379  h_Pkaon_den->Sumw2();
380  h_Pkaon_num->Sumw2();
381  h_Pmichel_e_den->Sumw2();
382  h_Pmichel_e_num->Sumw2();
383  h_Efrac_kaon = tfs->make<TH1D>("h_Efrac_kaon", "Efrac Kaon; Track Purity;", 60, 0, 1.2);
385  tfs->make<TH1D>("h_Ecomplet_kaon", "Ecomplet Kaon; Track Completeness;", 60, 0, 1.2);
386  h_trackRes_kaon = tfs->make<TH1D>(
387  "h_trackRes_kaon", "Kaon Residual; Truth length - Reco length (cm);", 200, -100, 100);
389  tfs->make<TH1D>("h_Efrac_michel", "Efrac Michel; Track Energy fraction;", 60, 0, 1.2);
391  tfs->make<TH1D>("h_Ecomplet_michel", "Ecomplet Michel; Track Completeness;", 60, 0, 1.2);
392  h_trackRes_michel = tfs->make<TH1D>(
393  "h_trackRes_michel", "Michel Residual; Truth length - Reco length (cm);", 200, -100, 100);
394  h_kaon_length =
395  tfs->make<TH1D>("h_kaon_length", "Kaon Length; Kaon Truth Length (cm)", 40, 0, 100);
397  tfs->make<TH1D>("h_kaonwtrk_length", "Kaon Length; Kaon Truth Length (cm)", 40, 0, 100);
399  tfs->make<TH1D>("h_michel_length", "Michel Length; Michel e Truth Length (cm)", 40, 0, 100);
400  h_michelwtrk_length = tfs->make<TH1D>(
401  "h_michelwtrk_length", "Michel Length; Michel e Truth Length (cm)", 40, 0, 100);
402 
403  h_Efrac_kaon->Sumw2();
404  h_Ecomplet_kaon->Sumw2();
405  h_trackRes_kaon->Sumw2();
406  h_Efrac_michel->Sumw2();
407  h_Ecomplet_michel->Sumw2();
408  h_trackRes_michel->Sumw2();
409  h_kaon_length->Sumw2();
410  h_kaonwtrk_length->Sumw2();
411  h_michel_length->Sumw2();
412  h_michelwtrk_length->Sumw2();
413  }
Geometry information for a single TPC.
Definition: TPCGeo.h:38
then local
art::ServiceHandle< art::TFileService > tfs
void LocalToWorld(const double *tpc, double *world) const
Transform point from local TPC frame to world frame.
Definition: TPCGeo.h:563
BEGIN_PROLOG could also be cout
void DUNE::NeutrinoTrackingEff::beginRun ( const art::Run &  run)
private

Definition at line 422 of file NeutrinoTrackingEff_module.cc.

423  {
424  mf::LogInfo("NeutrinoTrackingEff") << "begin run..." << std::endl;
425  }
void DUNE::NeutrinoTrackingEff::doEfficiencies ( )
private

Definition at line 1009 of file NeutrinoTrackingEff_module.cc.

1010  {
1011  art::ServiceHandle<art::TFileService const> tfs;
1012 
1013  if (TEfficiency::CheckConsistency(*h_Ev_num, *h_Ev_den)) {
1014  h_Eff_Ev = tfs->make<TEfficiency>(*h_Ev_num, *h_Ev_den);
1015  TGraphAsymmErrors* grEff_Ev = h_Eff_Ev->CreateGraph();
1016  grEff_Ev->Write("grEff_Ev");
1017  h_Eff_Ev->Write("h_Eff_Ev");
1018  }
1019  if (TEfficiency::CheckConsistency(*h_Pmu_num, *h_Pmu_den)) {
1020  h_Eff_Pmu = tfs->make<TEfficiency>(*h_Pmu_num, *h_Pmu_den);
1021  TGraphAsymmErrors* grEff_Pmu = h_Eff_Pmu->CreateGraph();
1022  grEff_Pmu->Write("grEff_Pmu");
1023  h_Eff_Pmu->Write("h_Eff_Pmu");
1024  }
1025  if (TEfficiency::CheckConsistency(*h_theta_num, *h_theta_den)) {
1026  h_Eff_theta = tfs->make<TEfficiency>(*h_theta_num, *h_theta_den);
1027  TGraphAsymmErrors* grEff_theta = h_Eff_theta->CreateGraph();
1028  grEff_theta->Write("grEff_theta");
1029  h_Eff_theta->Write("h_Eff_theta");
1030  }
1031  if (TEfficiency::CheckConsistency(*h_Pproton_num, *h_Pproton_den)) {
1032  h_Eff_Pproton = tfs->make<TEfficiency>(*h_Pproton_num, *h_Pproton_den);
1033  TGraphAsymmErrors* grEff_Pproton = h_Eff_Pproton->CreateGraph();
1034  grEff_Pproton->Write("grEff_Pproton");
1035  h_Eff_Pproton->Write("h_Eff_Pproton");
1036  }
1037  if (TEfficiency::CheckConsistency(*h_Ppion_plus_num, *h_Ppion_plus_den)) {
1038  h_Eff_Ppion_plus = tfs->make<TEfficiency>(*h_Ppion_plus_num, *h_Ppion_plus_den);
1039  TGraphAsymmErrors* grEff_Ppion_plus = h_Eff_Ppion_plus->CreateGraph();
1040  grEff_Ppion_plus->Write("grEff_Ppion_plus");
1041  h_Eff_Ppion_plus->Write("h_Eff_Ppion_plus");
1042  }
1043  if (TEfficiency::CheckConsistency(*h_Ppion_minus_num, *h_Ppion_minus_den)) {
1044  h_Eff_Ppion_minus = tfs->make<TEfficiency>(*h_Ppion_minus_num, *h_Ppion_minus_den);
1045  TGraphAsymmErrors* grEff_Ppion_minus = h_Eff_Ppion_minus->CreateGraph();
1046  grEff_Ppion_minus->Write("grEff_Ppion_minus");
1047  h_Eff_Ppion_minus->Write("h_Eff_Ppion_minus");
1048  }
1049  if (TEfficiency::CheckConsistency(*h_muonwtrk_length, *h_muon_length)) {
1050  h_Eff_Lmuon = tfs->make<TEfficiency>(*h_muonwtrk_length, *h_muon_length);
1051  TGraphAsymmErrors* grEff_Lmuon = h_Eff_Lmuon->CreateGraph();
1052  grEff_Lmuon->Write("grEff_Lmuon");
1053  h_Eff_Lmuon->Write("h_Eff_Lmuon");
1054  }
1055  if (TEfficiency::CheckConsistency(*h_protonwtrk_length, *h_proton_length)) {
1056  h_Eff_Lproton = tfs->make<TEfficiency>(*h_protonwtrk_length, *h_proton_length);
1057  TGraphAsymmErrors* grEff_Lproton = h_Eff_Lproton->CreateGraph();
1058  grEff_Lproton->Write("grEff_Lproton");
1059  h_Eff_Lproton->Write("h_Eff_Lproton");
1060  }
1061  if (TEfficiency::CheckConsistency(*h_pionpwtrk_length, *h_pionp_length)) {
1062  h_Eff_Lpion_plus = tfs->make<TEfficiency>(*h_pionpwtrk_length, *h_pionp_length);
1063  TGraphAsymmErrors* grEff_Lpion_plus = h_Eff_Lpion_plus->CreateGraph();
1064  grEff_Lpion_plus->Write("grEff_Lpion_plus");
1065  h_Eff_Lpion_plus->Write("h_Eff_Lpion_plus");
1066  }
1067  if (TEfficiency::CheckConsistency(*h_pionpwtrk_length, *h_pionp_length)) {
1068  h_Eff_Lpion_minus = tfs->make<TEfficiency>(*h_pionmwtrk_length, *h_pionm_length);
1069  TGraphAsymmErrors* grEff_Lpion_minus = h_Eff_Lpion_minus->CreateGraph();
1070  grEff_Lpion_minus->Write("grEff_Lpion_minus");
1071  h_Eff_Lpion_minus->Write("h_Eff_Lpion_minus");
1072  }
1073  if (TEfficiency::CheckConsistency(*h_Pkaon_num, *h_Pkaon_den)) {
1074  h_Eff_Pkaon = tfs->make<TEfficiency>(*h_Pkaon_num, *h_Pkaon_den);
1075  TGraphAsymmErrors* grEff_Pkaon = h_Eff_Pkaon->CreateGraph();
1076  grEff_Pkaon->Write("grEff_Pkaon");
1077  h_Eff_Pkaon->Write("h_Eff_Pkaon");
1078  }
1079  if (TEfficiency::CheckConsistency(*h_kaonwtrk_length, *h_kaon_length)) {
1080  h_Eff_Lkaon = tfs->make<TEfficiency>(*h_kaonwtrk_length, *h_kaon_length);
1081  TGraphAsymmErrors* grEff_Lkaon = h_Eff_Lkaon->CreateGraph();
1082  grEff_Lkaon->Write("grEff_Lkaon");
1083  h_Eff_Lkaon->Write("h_Eff_Lkaon");
1084  }
1085  if (TEfficiency::CheckConsistency(*h_Pmichel_e_num, *h_Pmichel_e_den)) {
1086  h_Eff_Pmichel = tfs->make<TEfficiency>(*h_Pmichel_e_num, *h_Pmichel_e_den);
1087  TGraphAsymmErrors* grEff_Pmichel = h_Eff_Pmichel->CreateGraph();
1088  grEff_Pmichel->Write("grEff_Pmichel");
1089  h_Eff_Pmichel->Write("h_Eff_Pmichel");
1090  }
1091  if (TEfficiency::CheckConsistency(*h_michelwtrk_length, *h_michel_length)) {
1092  h_Eff_Lmichel = tfs->make<TEfficiency>(*h_michelwtrk_length, *h_michel_length);
1093  TGraphAsymmErrors* grEff_Lmichel = h_Eff_Lmichel->CreateGraph();
1094  grEff_Lmichel->Write("grEff_Lmichel");
1095  h_Eff_Lmichel->Write("h_Eff_Lmichel");
1096  }
1097  }
art::ServiceHandle< art::TFileService > tfs
void DUNE::NeutrinoTrackingEff::endJob ( )
private

Definition at line 416 of file NeutrinoTrackingEff_module.cc.

417  {
418  doEfficiencies();
419  }
bool DUNE::NeutrinoTrackingEff::insideFV ( double  vertex[4])
private

Definition at line 998 of file NeutrinoTrackingEff_module.cc.

999  {
1000  double const x = vertex[0];
1001  double const y = vertex[1];
1002  double const z = vertex[2];
1003 
1004  return x > fFidVolXmin && x < fFidVolXmax && y > fFidVolYmin && y < fFidVolYmax &&
1005  z > fFidVolZmin && z < fFidVolZmax;
1006  }
process_name vertex
Definition: cheaterreco.fcl:51
process_name opflash particleana ie ie ie z
process_name opflash particleana ie x
process_name opflash particleana ie ie y
void DUNE::NeutrinoTrackingEff::processEff ( const art::Event &  evt)
private

save FS particles

Definition at line 436 of file NeutrinoTrackingEff_module.cc.

437  {
438  // Save neutrino's interaction info
439  art::Handle<std::vector<simb::MCTruth>> MCtruthHandle;
440  event.getByLabel(fMCTruthModuleLabel, MCtruthHandle);
441  std::vector<art::Ptr<simb::MCTruth>> MCtruthlist;
442  art::fill_ptr_vector(MCtruthlist, MCtruthHandle);
443  art::Ptr<simb::MCTruth> MCtruth;
444 
445  //For now assume that there is only one neutrino interaction...
446  MCtruth = MCtruthlist[0];
447  if (MCtruth->NeutrinoSet()) {
448  simb::MCNeutrino nu = MCtruth->GetNeutrino();
449  if (nu.CCNC() == 0)
450  MC_isCC = 1;
451  else if (nu.CCNC() == 1)
452  MC_isCC = 0;
453  simb::MCParticle neutrino = nu.Nu();
454  MC_incoming_PDG = nu.Nu().PdgCode();
455  const TLorentzVector& nu_momentum = nu.Nu().Momentum(0);
456  nu_momentum.GetXYZT(MC_incoming_P);
457  const TLorentzVector& vertex = neutrino.Position(0);
458  vertex.GetXYZT(MC_vertex);
459  }
460 
461  //!save FS particles
462 
463  double tmp_leadingPionPlusE = 0.0;
464  double tmp_leadingPionMinusE = 0.0;
465  double tmp_leadingProtonE = 0.0;
466 
467  simb::MCParticle* MClepton = nullptr;
468  simb::MCParticle* MCproton = nullptr;
469  simb::MCParticle* MCpion_plus = nullptr;
470  simb::MCParticle* MCpion_minus = nullptr;
471  simb::MCParticle* MCkaon = nullptr;
472  simb::MCParticle* MCmichel = nullptr;
473 
474  art::ServiceHandle<cheat::ParticleInventoryService const> pi_serv;
475  const sim::ParticleList& plist = pi_serv->ParticleList();
476  simb::MCParticle* particle = 0;
477  int i = 0; // particle index
478 
479  auto const clockData =
480  art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(event);
481  auto const detProp =
482  art::ServiceHandle<detinfo::DetectorPropertiesService const>()->DataFor(event);
483 
484  for (sim::ParticleList::const_iterator ipar = plist.begin(); ipar != plist.end(); ++ipar) {
485  particle = ipar->second;
486  if (particle->PdgCode() == fLeptonPDGcode && particle->Mother() == 0) { //primary lepton
487  const TLorentzVector& lepton_momentum = particle->Momentum(0);
488  lepton_momentum.GetXYZT(MC_lepton_startMomentum);
489  MC_leptonID = particle->TrackId();
490  MC_leptonP = sqrt(pow(particle->Momentum().Px(), 2) + pow(particle->Momentum().Py(), 2) +
491  pow(particle->Momentum().Pz(), 2));
492  MClepton = particle;
493  }
494  if (particle->Mother() == 0) { //save primary particle i.e. from the neutrino interaction
495  //save leading pion and proton
496  if (particle->PdgCode() == 2212) {
497  if (particle->Momentum().E() > tmp_leadingProtonE) {
498  tmp_leadingProtonE = particle->Momentum().E();
499  MC_leading_protonID = particle->TrackId();
501  sqrt(pow(particle->Momentum().Px(), 2) + pow(particle->Momentum().Py(), 2) +
502  pow(particle->Momentum().Pz(), 2));
503  MCproton = particle;
504  }
505  }
506  else if (particle->PdgCode() == 211) {
507  if (particle->Momentum().E() > tmp_leadingPionPlusE) {
508  tmp_leadingPionPlusE = particle->Momentum().E();
509  MC_leading_PionPlusID = particle->TrackId();
511  sqrt(pow(particle->Momentum().Px(), 2) + pow(particle->Momentum().Py(), 2) +
512  pow(particle->Momentum().Pz(), 2));
513  MCpion_plus = particle;
514  }
515  }
516  else if (particle->PdgCode() == -211) {
517  if (particle->Momentum().E() > tmp_leadingPionMinusE) {
518  tmp_leadingPionMinusE = particle->Momentum().E();
519  MC_leading_PionMinusID = particle->TrackId();
521  sqrt(pow(particle->Momentum().Px(), 2) + pow(particle->Momentum().Py(), 2) +
522  pow(particle->Momentum().Pz(), 2));
523  MCpion_minus = particle;
524  }
525  }
526  i++; //paticle index
527  }
528 
529  //=======================================================================================
530  //add Nucleon decay stuff and particle cannon
531  //=======================================================================================
532  if (!fisNeutrinoInt) {
533  if (particle->Mother() == 0) {
534  const TLorentzVector& positionStart = particle->Position(0);
535  positionStart.GetXYZT(MC_vertex);
536  }
537  if (particle->PdgCode() == 321) { //save primary Kaon
538  MC_kaonID = particle->TrackId();
539  MC_kaonP = sqrt(pow(particle->Momentum().Px(), 2) + pow(particle->Momentum().Py(), 2) +
540  pow(particle->Momentum().Pz(), 2));
541  MCkaon = particle;
542  }
543  else if (particle->PdgCode() == fLeptonPDGcode) { // Particle cannon muon
544  MC_leptonID = particle->TrackId();
545  MC_leptonP = sqrt(pow(particle->Momentum().Px(), 2) + pow(particle->Momentum().Py(), 2) +
546  pow(particle->Momentum().Pz(), 2));
547  MClepton = particle;
548  }
549  else if (particle->PdgCode() == 2212) {
550  if (particle->Momentum().E() > tmp_leadingProtonE) {
551  tmp_leadingProtonE = particle->Momentum().E();
552  MC_leading_protonID = particle->TrackId();
554  sqrt(pow(particle->Momentum().Px(), 2) + pow(particle->Momentum().Py(), 2) +
555  pow(particle->Momentum().Pz(), 2));
556  MCproton = particle;
557  }
558  }
559  else if (particle->PdgCode() == 211) {
560  if (particle->Momentum().E() > tmp_leadingPionPlusE) {
561  tmp_leadingPionPlusE = particle->Momentum().E();
562  MC_leading_PionPlusID = particle->TrackId();
564  sqrt(pow(particle->Momentum().Px(), 2) + pow(particle->Momentum().Py(), 2) +
565  pow(particle->Momentum().Pz(), 2));
566  MCpion_plus = particle;
567  }
568  }
569  else if (particle->PdgCode() == -211) {
570  if (particle->Momentum().E() > tmp_leadingPionMinusE) {
571  tmp_leadingPionMinusE = particle->Momentum().E();
572  MC_leading_PionMinusID = particle->TrackId();
574  sqrt(pow(particle->Momentum().Px(), 2) + pow(particle->Momentum().Py(), 2) +
575  pow(particle->Momentum().Pz(), 2));
576  MCpion_minus = particle;
577  }
578  }
579  else if (particle->Process() == "Decay" &&
580  particle->PdgCode() == -11) { // michel electron from muon decay
581  MC_michelID = particle->TrackId();
582  MC_michelP = sqrt(pow(particle->Momentum().Px(), 2) + pow(particle->Momentum().Py(), 2) +
583  pow(particle->Momentum().Pz(), 2));
584  MCmichel = particle;
585  }
586  else if (TMath::Abs(particle->PdgCode() == 321)) { //save primary Kaon
587  MC_kaonID = particle->TrackId();
588  MC_kaonP = sqrt(pow(particle->Momentum().Px(), 2) + pow(particle->Momentum().Py(), 2) +
589  pow(particle->Momentum().Pz(), 2));
590  MCkaon = particle;
591  }
592  }
593  }
594  //===================================================================
595  //Saving denominator histograms
596  //===================================================================
597  if (not insideFV(MC_vertex)) return;
598  double Pv =
599  sqrt(pow(MC_incoming_P[0], 2) + pow(MC_incoming_P[1], 2) + pow(MC_incoming_P[2], 2));
600  double theta_mu = acos((MC_incoming_P[0] * MC_lepton_startMomentum[0] +
603  (Pv * MC_leptonP));
604  theta_mu *= (180.0 / 3.14159);
605  double truth_lengthLepton = truthLength(clockData, detProp, MClepton);
606  double proton_length = truthLength(clockData, detProp, MCproton);
607  double pion_plus_length = truthLength(clockData, detProp, MCpion_plus);
608  double pion_minus_length = truthLength(clockData, detProp, MCpion_minus);
609  double kaonLength = truthLength(clockData, detProp, MCkaon);
610  double michelLength = truthLength(clockData, detProp, MCmichel);
611 
612  // save CC events within the fiducial volume with the favorite neutrino
613  // flavor
615  if (MClepton) {
616  h_Ev_den->Fill(MC_incoming_P[3]);
617  h_Pmu_den->Fill(MC_leptonP);
618  h_theta_den->Fill(theta_mu);
619  h_muon_length->Fill(truth_lengthLepton);
620  }
621  if (MCproton) {
623  h_proton_length->Fill(proton_length);
624  }
625  if (MCpion_plus) {
627  h_pionp_length->Fill(pion_plus_length);
628  }
629  if (MCpion_minus) {
631  h_pionm_length->Fill(pion_minus_length);
632  }
633  if (MCkaon) {
634  h_Pkaon_den->Fill(MC_kaonP);
635  h_kaon_length->Fill(kaonLength);
636  }
637  }
638 
639  //save events for Nucleon decay and particle cannon
640  if (!fisNeutrinoInt) {
641  if (MClepton) {
642  h_Pmu_den->Fill(MC_leptonP);
643  h_muon_length->Fill(truth_lengthLepton);
644  }
645  if (MCkaon) {
646  h_Pkaon_den->Fill(MC_kaonP);
647  h_kaon_length->Fill(kaonLength);
648  }
649  if (MCproton) {
651  h_proton_length->Fill(proton_length);
652  }
653  if (MCpion_plus) {
655  h_pionp_length->Fill(pion_plus_length);
656  }
657  if (MCpion_minus) {
659  h_pionm_length->Fill(pion_minus_length);
660  }
661  if (MCmichel) {
663  h_michel_length->Fill(michelLength);
664  }
665  }
666 
667  //========================================================================
668  // Reco stuff, once we have selected a MC Particle let's find out if there is a track associated
669  //========================================================================
670  art::Handle<std::vector<recob::Track>> trackListHandle;
671  if (!event.getByLabel(fTrackModuleLabel, trackListHandle)) return;
672  std::vector<art::Ptr<recob::Track>> tracklist;
673  art::fill_ptr_vector(tracklist, trackListHandle);
674  int n_recoTrack = tracklist.size();
675 
676  art::FindManyP<recob::Hit> track_hits(trackListHandle, event, fTrackModuleLabel);
677  if (n_recoTrack == 0) {
678  MF_LOG_DEBUG("NeutrinoTrackingEff") << "There are no reco tracks... bye";
679  return;
680  }
681  MF_LOG_DEBUG("NeutrinoTrackingEff") << "Found this many reco tracks " << n_recoTrack;
682 
683  double Efrac_lepton = 0.0;
684  double Ecomplet_lepton = 0.0;
685  double Efrac_proton = 0.0;
686  double Ecomplet_proton = 0.0;
687  double Efrac_pionplus = 0.0;
688  double Ecomplet_pionplus = 0.0;
689  double Efrac_pionminus = 0.0;
690  double Ecomplet_pionminus = 0.0;
691  double Efrac_kaon = 0.0;
692  double Ecomplet_kaon = 0.0;
693  double Efrac_michel = 0.0;
694  double Ecomplet_michel = 0.0;
695  double trackLength_lepton = 0.0;
696  double trackLength_proton = 0.0;
697  double trackLength_pion_plus = 0.0;
698  double trackLength_pion_minus = 0.0;
699  double trackLength_kaon = 0.0;
700  double trackLength_michel = 0.0;
701  const simb::MCParticle* MClepton_reco = nullptr;
702  const simb::MCParticle* MCproton_reco = nullptr;
703  const simb::MCParticle* MCpion_plus_reco = nullptr;
704  const simb::MCParticle* MCpion_minus_reco = nullptr;
705  const simb::MCParticle* MCkaon_reco = nullptr;
706  const simb::MCParticle* MCmichel_reco = nullptr;
707 
708  std::vector<art::Ptr<recob::Hit>> tmp_all_trackHits = track_hits.at(0);
709  std::vector<art::Ptr<recob::Hit>> all_hits;
710  art::Handle<std::vector<recob::Hit>> hithandle;
711  auto const pd = event.getProductDescription(tmp_all_trackHits[0].id());
712  if (pd && event.getByLabel(pd->inputTag(), hithandle)) {
713  art::fill_ptr_vector(all_hits, hithandle);
714  }
715 
716  for (int i = 0; i < n_recoTrack; i++) {
717  art::Ptr<recob::Track> track = tracklist[i];
718  std::vector<art::Ptr<recob::Hit>> all_trackHits = track_hits.at(i);
719  double tmpEfrac = 0;
720  double tmpEcomplet = 0;
721  const simb::MCParticle* particle;
722  truthMatcher(clockData, all_hits, all_trackHits, particle, tmpEfrac, tmpEcomplet);
723  if (!particle) continue;
724  if ((particle->PdgCode() == fLeptonPDGcode) && (particle->TrackId() == MC_leptonID)) {
725  // save the best track ... based on completeness if there is more than
726  // one track if( tmpEfrac > Efrac_lepton ){ ///this was base on purity
727  if (tmpEcomplet > Ecomplet_lepton) {
728  Ecomplet_lepton = tmpEcomplet;
729  Efrac_lepton = tmpEfrac;
730  trackLength_lepton = track->Length();
731  MClepton_reco = particle;
732  }
733  }
734  else if ((particle->PdgCode() == 2212) && (particle->TrackId() == MC_leading_protonID)) {
735  //save the best track ... based on completeness if there is more than one track
736  if (tmpEcomplet > Ecomplet_proton) {
737  Ecomplet_proton = tmpEcomplet;
738  Efrac_proton = tmpEfrac;
739  trackLength_proton = track->Length();
740  MCproton_reco = particle;
741  }
742  }
743  else if ((particle->PdgCode() == 211) && (particle->TrackId() == MC_leading_PionPlusID)) {
744  //save the best track ... based on completeness if there is more than one track
745  if (tmpEcomplet > Ecomplet_pionplus) {
746  Ecomplet_pionplus = tmpEcomplet;
747  Efrac_pionplus = tmpEfrac;
748  trackLength_pion_plus = track->Length();
749  MCpion_plus_reco = particle;
750  }
751  }
752  else if ((particle->PdgCode() == -211) && (particle->TrackId() == MC_leading_PionMinusID)) {
753  //save the best track ... based on completeness if there is more than one track
754  if (tmpEcomplet > Ecomplet_pionminus) {
755  Ecomplet_pionminus = tmpEcomplet;
756  Efrac_pionminus = tmpEfrac;
757  trackLength_pion_minus = track->Length();
758  MCpion_minus_reco = particle;
759  }
760  }
761  //kaon from nucleon decay
762  else if ((TMath::Abs(particle->PdgCode()) == 321) && (particle->TrackId() == MC_kaonID)) {
763  //save the best track ... based on completeness if there is more than one track
764  if (tmpEcomplet > Ecomplet_kaon) {
765  Ecomplet_kaon = tmpEcomplet;
766  Efrac_kaon = tmpEfrac;
767  trackLength_kaon = track->Length();
768  MCkaon_reco = particle;
769  }
770  }
771  //michel from nucleon decay
772  else if ((particle->PdgCode() == -11) && (particle->TrackId() == MC_michelID)) {
773  //save the best track ... based on completeness if there is more than one track
774  if (tmpEcomplet > Ecomplet_michel) {
775  Ecomplet_michel = tmpEcomplet;
776  Efrac_michel = tmpEfrac;
777  trackLength_michel = track->Length();
778  MCmichel_reco = particle;
779  }
780  }
781  }
782 
783  double Reco_LengthRes = truth_lengthLepton - trackLength_lepton;
784  double Reco_LengthResProton = proton_length - trackLength_proton;
785  double Reco_LengthResPionPlus = pion_plus_length - trackLength_pion_plus;
786  double Reco_LengthResPionMinus = pion_minus_length - trackLength_pion_minus;
787 
788  if (MClepton_reco && MClepton) {
790  h_Pmu_num->Fill(MC_leptonP);
791  h_Ev_num->Fill(MC_incoming_P[3]);
792  h_theta_num->Fill(theta_mu);
793  h_Efrac_lepton->Fill(Efrac_lepton);
794  h_Ecomplet_lepton->Fill(Ecomplet_lepton);
795  h_trackRes_lepton->Fill(Reco_LengthRes);
796  h_muonwtrk_length->Fill(truth_lengthLepton);
797  }
798  }
799  if (MCproton_reco && MCproton) {
802  h_Efrac_proton->Fill(Efrac_proton);
803  h_Ecomplet_proton->Fill(Ecomplet_proton);
804  h_trackRes_proton->Fill(Reco_LengthResProton);
805  h_protonwtrk_length->Fill(proton_length);
806  }
807  }
808  if (MCpion_plus_reco && MCpion_plus) {
811  h_Efrac_pion_plus->Fill(Efrac_pionplus);
812  h_Ecomplet_pion_plus->Fill(Ecomplet_pionplus);
813  h_trackRes_pion_plus->Fill(Reco_LengthResPionPlus);
814  h_pionpwtrk_length->Fill(pion_plus_length);
815  }
816  }
817  if (MCpion_minus_reco && MCpion_minus) {
820  h_Efrac_pion_minus->Fill(Efrac_pionminus);
821  h_Ecomplet_pion_minus->Fill(Ecomplet_pionminus);
822  h_trackRes_pion_minus->Fill(Reco_LengthResPionMinus);
823  h_pionmwtrk_length->Fill(pion_minus_length);
824  }
825  }
826  if (MCkaon_reco && MCkaon) {
828  h_Pkaon_num->Fill(MC_kaonP);
829  h_Efrac_kaon->Fill(Efrac_kaon);
830  h_Ecomplet_kaon->Fill(Ecomplet_kaon);
831  h_trackRes_kaon->Fill(kaonLength - trackLength_kaon);
832  h_kaonwtrk_length->Fill(kaonLength);
833  }
834  }
835  //Non neutrino events
836  //=========================================================
837  if (!fisNeutrinoInt) {
838  if (MClepton_reco && MClepton) {
839  h_Pmu_num->Fill(MC_leptonP);
840  h_Efrac_lepton->Fill(Efrac_lepton);
841  h_Ecomplet_lepton->Fill(Ecomplet_lepton);
842  h_trackRes_lepton->Fill(Reco_LengthRes);
843  h_muonwtrk_length->Fill(truth_lengthLepton);
844  }
845  if (MCkaon_reco && MCkaon) {
846  h_Pkaon_num->Fill(MC_kaonP);
847  h_Efrac_kaon->Fill(Efrac_kaon);
848  h_Ecomplet_kaon->Fill(Ecomplet_kaon);
849  h_trackRes_kaon->Fill(kaonLength - trackLength_kaon);
850  h_kaonwtrk_length->Fill(kaonLength);
851  }
852  if (MCproton_reco && MCproton) {
854  h_Efrac_proton->Fill(Efrac_proton);
855  h_Ecomplet_proton->Fill(Ecomplet_proton);
856  h_trackRes_proton->Fill(Reco_LengthResProton);
857  h_protonwtrk_length->Fill(proton_length);
858  }
859  if (MCpion_plus_reco && MCpion_plus) {
861  h_Efrac_pion_plus->Fill(Efrac_pionplus);
862  h_Ecomplet_pion_plus->Fill(Ecomplet_pionplus);
863  h_trackRes_pion_plus->Fill(Reco_LengthResPionPlus);
864  h_pionpwtrk_length->Fill(pion_plus_length);
865  }
866  if (MCpion_minus_reco && MCpion_minus) {
868  h_Efrac_pion_minus->Fill(Efrac_pionminus);
869  h_Ecomplet_pion_minus->Fill(Ecomplet_pionminus);
870  h_trackRes_pion_minus->Fill(Reco_LengthResPionMinus);
871  h_pionmwtrk_length->Fill(pion_minus_length);
872  }
873  if (MCmichel_reco && MCmichel) {
875  h_Efrac_michel->Fill(Efrac_michel);
876  h_Ecomplet_michel->Fill(Ecomplet_michel);
877  h_trackRes_michel->Fill(michelLength - trackLength_michel);
878  h_michelwtrk_length->Fill(michelLength);
879  }
880  }
881  }
process_name vertex
Definition: cheaterreco.fcl:51
process_name use argoneut_mc_hitfinder track
void truthMatcher(detinfo::DetectorClocksData const &clockData, std::vector< art::Ptr< recob::Hit >> all_hits, std::vector< art::Ptr< recob::Hit >> track_hits, const simb::MCParticle *&MCparticle, double &Efrac, double &Ecomplet)
double truthLength(const detinfo::DetectorClocksData &clockData, detinfo::DetectorPropertiesData const &detProp, const simb::MCParticle *MCparticle)
BEGIN_PROLOG SN nu
auto const detProp
double DUNE::NeutrinoTrackingEff::truthLength ( const detinfo::DetectorClocksData clockData,
detinfo::DetectorPropertiesData const &  detProp,
const simb::MCParticle *  MCparticle 
)
private

Definition at line 949 of file NeutrinoTrackingEff_module.cc.

952  {
953  // Calculate the truth length considering only the part that is
954  // inside the TPC Base on a peace of code from
955  // dune/TrackingAna/TrackingEfficiency_module.cc
956 
957  if (!MCparticle) return -999.0;
958  int numberTrajectoryPoints = MCparticle->NumberTrajectoryPoints();
959  std::vector<double> TPCLengthHits(numberTrajectoryPoints, 0.0);
960  int FirstHit = 0, LastHit = 0;
961  double TPCLength = 0.0;
962  bool BeenInVolume = false;
963 
964  double const WindowSize = detProp.NumberTimeSamples() * clockData.TPCClock().TickPeriod() * 1e3;
965 
966  for (unsigned int MCHit = 0; MCHit < TPCLengthHits.size(); ++MCHit) {
967  const TLorentzVector& tmpPosition = MCparticle->Position(MCHit);
968  double const tmpPosArray[] = {tmpPosition[0], tmpPosition[1], tmpPosition[2]};
969  if (MCHit != 0)
970  TPCLengthHits[MCHit] = sqrt(pow((MCparticle->Vx(MCHit - 1) - MCparticle->Vx(MCHit)), 2) +
971  pow((MCparticle->Vy(MCHit - 1) - MCparticle->Vy(MCHit)), 2) +
972  pow((MCparticle->Vz(MCHit - 1) - MCparticle->Vz(MCHit)), 2));
973  geo::TPCID tpcid = geom->FindTPCAtPosition(tmpPosArray);
974  if (tpcid.isValid) {
975  // -- Check if hit is within drift window...
976  geo::CryostatGeo const& cryo = geom->Cryostat(tpcid.Cryostat);
977  geo::TPCGeo const& tpc = cryo.TPC(tpcid.TPC);
978  double XPlanePosition = tpc.PlaneLocation(0)[0];
979  double DriftTimeCorrection = fabs(tmpPosition[0] - XPlanePosition) / fDriftVelocity;
980  double TimeAtPlane = MCparticle->T() + DriftTimeCorrection;
981 
982  if (TimeAtPlane < trigger_offset(clockData) ||
983  TimeAtPlane > trigger_offset(clockData) + WindowSize)
984  continue;
985  LastHit = MCHit;
986  if (!BeenInVolume) {
987  BeenInVolume = true;
988  FirstHit = MCHit;
989  }
990  }
991  }
992  for (int Hit = FirstHit + 1; Hit <= LastHit; ++Hit)
993  TPCLength += TPCLengthHits[Hit];
994  return TPCLength;
995  }
art::ServiceHandle< geo::Geometry const > geom
ElecClock const & TPCClock() const noexcept
Borrow a const TPC clock with time set to Trigger time [us].
Geometry information for a single TPC.
Definition: TPCGeo.h:38
Geometry information for a single cryostat.
Definition: CryostatGeo.h:43
constexpr double TickPeriod() const noexcept
A single tick period in microseconds.
Definition: ElecClock.h:352
The data type to uniquely identify a TPC.
Definition: geo_types.h:386
BEGIN_PROLOG triggeremu_data_config_icarus settings sequence::triggeremu_data_config_icarus settings PMTADCthresholds sequence::triggeremu_data_config_icarus settings PMTADCthresholds WindowSize
const TPCGeo & TPC(unsigned int itpc) const
Return the itpc&#39;th TPC in the cryostat.
Definition: CryostatGeo.cxx:93
int trigger_offset(DetectorClocksData const &data)
const double * PlaneLocation(unsigned int p) const
Definition: TPCGeo.cxx:382
auto const detProp
void DUNE::NeutrinoTrackingEff::truthMatcher ( detinfo::DetectorClocksData const &  clockData,
std::vector< art::Ptr< recob::Hit >>  all_hits,
std::vector< art::Ptr< recob::Hit >>  track_hits,
const simb::MCParticle *&  MCparticle,
double &  Efrac,
double &  Ecomplet 
)
private

Definition at line 884 of file NeutrinoTrackingEff_module.cc.

890  {
891  art::ServiceHandle<cheat::BackTrackerService const> bt_serv;
892  art::ServiceHandle<cheat::ParticleInventoryService const> pi_serv;
893  std::map<int, double> trkID_E;
894  for (size_t j = 0; j < track_hits.size(); ++j) {
895  art::Ptr<recob::Hit> hit = track_hits[j];
896  std::vector<sim::TrackIDE> TrackIDs = bt_serv->HitToTrackIDEs(clockData, hit);
897  for (size_t k = 0; k < TrackIDs.size(); k++) {
898  trkID_E[TrackIDs[k].trackID] += TrackIDs[k].energy;
899  }
900  }
901  double E_em = 0.0;
902  double max_E = -999.0;
903  double total_E = 0.0;
904  int TrackID = -999;
905  double partial_E =
906  0.0; // amount of energy deposited by the particle that deposited more energy...
907 
908  // If the collection of hits have more than one particle associate
909  // save the particle w/ the highest energy deposition since we are
910  // looking for muons/pions/protons this should be enough
911  if (!trkID_E.size()) {
912  MCparticle = 0;
913  return; //Ghost track???
914  }
915  for (std::map<int, double>::iterator ii = trkID_E.begin(); ii != trkID_E.end(); ++ii) {
916  total_E += ii->second;
917  if ((ii->second) > max_E) {
918  partial_E = ii->second;
919  max_E = ii->second;
920  TrackID = ii->first;
921  if (TrackID < 0) E_em += ii->second;
922  }
923  }
924  MCparticle = pi_serv->TrackIdToParticle_P(TrackID);
925 
926  // In the current simulation, we do not save EM Shower daughters
927  // in GEANT. But we do save the energy deposition in TrackIDEs. If
928  // the energy deposition is from a particle that is the daughter
929  // of an EM particle, the negative of the parent track ID is saved
930  // in TrackIDE for the daughter particle we don't want to track
931  // gammas or any other EM activity
932  if (TrackID < 0) return;
933 
934  Efrac = (partial_E) / total_E;
935 
936  // Completeness
937  double totenergy = 0;
938  for (size_t k = 0; k < all_hits.size(); ++k) {
939  art::Ptr<recob::Hit> hit = all_hits[k];
940  std::vector<sim::TrackIDE> TrackIDs = bt_serv->HitToTrackIDEs(clockData, hit);
941  for (size_t l = 0; l < TrackIDs.size(); ++l) {
942  if (TrackIDs[l].trackID == TrackID) totenergy += TrackIDs[l].energy;
943  }
944  }
945  Ecomplet = partial_E / totenergy;
946  }
std::vector< TrackID > TrackIDs
process_name hit
Definition: cheaterreco.fcl:51
pdgs k
Definition: selectors.fcl:22

Member Data Documentation

double DUNE::NeutrinoTrackingEff::fDriftVelocity
private

Definition at line 181 of file NeutrinoTrackingEff_module.cc.

float DUNE::NeutrinoTrackingEff::fFidVolCutX
private

Definition at line 170 of file NeutrinoTrackingEff_module.cc.

float DUNE::NeutrinoTrackingEff::fFidVolCutY
private

Definition at line 171 of file NeutrinoTrackingEff_module.cc.

float DUNE::NeutrinoTrackingEff::fFidVolCutZ
private

Definition at line 172 of file NeutrinoTrackingEff_module.cc.

float DUNE::NeutrinoTrackingEff::fFidVolXmax
private

Definition at line 175 of file NeutrinoTrackingEff_module.cc.

float DUNE::NeutrinoTrackingEff::fFidVolXmin
private

Definition at line 174 of file NeutrinoTrackingEff_module.cc.

float DUNE::NeutrinoTrackingEff::fFidVolYmax
private

Definition at line 177 of file NeutrinoTrackingEff_module.cc.

float DUNE::NeutrinoTrackingEff::fFidVolYmin
private

Definition at line 176 of file NeutrinoTrackingEff_module.cc.

float DUNE::NeutrinoTrackingEff::fFidVolZmax
private

Definition at line 179 of file NeutrinoTrackingEff_module.cc.

float DUNE::NeutrinoTrackingEff::fFidVolZmin
private

Definition at line 178 of file NeutrinoTrackingEff_module.cc.

bool DUNE::NeutrinoTrackingEff::fisNeutrinoInt
private

Definition at line 81 of file NeutrinoTrackingEff_module.cc.

int DUNE::NeutrinoTrackingEff::fLeptonPDGcode
private

Definition at line 78 of file NeutrinoTrackingEff_module.cc.

double DUNE::NeutrinoTrackingEff::fMaxLeptonP
private

Definition at line 80 of file NeutrinoTrackingEff_module.cc.

double DUNE::NeutrinoTrackingEff::fMaxNeutrinoE
private

Definition at line 79 of file NeutrinoTrackingEff_module.cc.

std::string DUNE::NeutrinoTrackingEff::fMCTruthModuleLabel
private

Definition at line 75 of file NeutrinoTrackingEff_module.cc.

int DUNE::NeutrinoTrackingEff::fNeutrinoPDGcode
private

Definition at line 77 of file NeutrinoTrackingEff_module.cc.

std::string DUNE::NeutrinoTrackingEff::fTrackModuleLabel
private

Definition at line 76 of file NeutrinoTrackingEff_module.cc.

art::ServiceHandle<geo::Geometry const> DUNE::NeutrinoTrackingEff::geom
private

Definition at line 182 of file NeutrinoTrackingEff_module.cc.

TH1D* DUNE::NeutrinoTrackingEff::h_Ecomplet_kaon
private

Definition at line 156 of file NeutrinoTrackingEff_module.cc.

TH1D* DUNE::NeutrinoTrackingEff::h_Ecomplet_lepton
private

Definition at line 117 of file NeutrinoTrackingEff_module.cc.

TH1D* DUNE::NeutrinoTrackingEff::h_Ecomplet_michel
private

Definition at line 159 of file NeutrinoTrackingEff_module.cc.

TH1D* DUNE::NeutrinoTrackingEff::h_Ecomplet_pion_minus
private

Definition at line 123 of file NeutrinoTrackingEff_module.cc.

TH1D* DUNE::NeutrinoTrackingEff::h_Ecomplet_pion_plus
private

Definition at line 121 of file NeutrinoTrackingEff_module.cc.

TH1D* DUNE::NeutrinoTrackingEff::h_Ecomplet_proton
private

Definition at line 119 of file NeutrinoTrackingEff_module.cc.

TEfficiency* DUNE::NeutrinoTrackingEff::h_Eff_Ev = 0
private

Definition at line 138 of file NeutrinoTrackingEff_module.cc.

TEfficiency* DUNE::NeutrinoTrackingEff::h_Eff_Lkaon = 0
private

Definition at line 167 of file NeutrinoTrackingEff_module.cc.

TEfficiency* DUNE::NeutrinoTrackingEff::h_Eff_Lmichel = 0
private

Definition at line 168 of file NeutrinoTrackingEff_module.cc.

TEfficiency* DUNE::NeutrinoTrackingEff::h_Eff_Lmuon = 0
private

Definition at line 145 of file NeutrinoTrackingEff_module.cc.

TEfficiency* DUNE::NeutrinoTrackingEff::h_Eff_Lpion_minus = 0
private

Definition at line 148 of file NeutrinoTrackingEff_module.cc.

TEfficiency* DUNE::NeutrinoTrackingEff::h_Eff_Lpion_plus = 0
private

Definition at line 147 of file NeutrinoTrackingEff_module.cc.

TEfficiency* DUNE::NeutrinoTrackingEff::h_Eff_Lproton = 0
private

Definition at line 146 of file NeutrinoTrackingEff_module.cc.

TEfficiency* DUNE::NeutrinoTrackingEff::h_Eff_Pkaon = 0
private

Definition at line 165 of file NeutrinoTrackingEff_module.cc.

TEfficiency* DUNE::NeutrinoTrackingEff::h_Eff_Pmichel = 0
private

Definition at line 166 of file NeutrinoTrackingEff_module.cc.

TEfficiency* DUNE::NeutrinoTrackingEff::h_Eff_Pmu = 0
private

Definition at line 139 of file NeutrinoTrackingEff_module.cc.

TEfficiency* DUNE::NeutrinoTrackingEff::h_Eff_Ppion_minus = 0
private

Definition at line 143 of file NeutrinoTrackingEff_module.cc.

TEfficiency* DUNE::NeutrinoTrackingEff::h_Eff_Ppion_plus = 0
private

Definition at line 142 of file NeutrinoTrackingEff_module.cc.

TEfficiency* DUNE::NeutrinoTrackingEff::h_Eff_Pproton = 0
private

Definition at line 141 of file NeutrinoTrackingEff_module.cc.

TEfficiency* DUNE::NeutrinoTrackingEff::h_Eff_theta = 0
private

Definition at line 140 of file NeutrinoTrackingEff_module.cc.

TH1D* DUNE::NeutrinoTrackingEff::h_Efrac_kaon
private

Definition at line 155 of file NeutrinoTrackingEff_module.cc.

TH1D* DUNE::NeutrinoTrackingEff::h_Efrac_lepton
private

Definition at line 116 of file NeutrinoTrackingEff_module.cc.

TH1D* DUNE::NeutrinoTrackingEff::h_Efrac_michel
private

Definition at line 158 of file NeutrinoTrackingEff_module.cc.

TH1D* DUNE::NeutrinoTrackingEff::h_Efrac_pion_minus
private

Definition at line 122 of file NeutrinoTrackingEff_module.cc.

TH1D* DUNE::NeutrinoTrackingEff::h_Efrac_pion_plus
private

Definition at line 120 of file NeutrinoTrackingEff_module.cc.

TH1D* DUNE::NeutrinoTrackingEff::h_Efrac_proton
private

Definition at line 118 of file NeutrinoTrackingEff_module.cc.

TH1D* DUNE::NeutrinoTrackingEff::h_Ev_den
private

Definition at line 103 of file NeutrinoTrackingEff_module.cc.

TH1D* DUNE::NeutrinoTrackingEff::h_Ev_num
private

Definition at line 104 of file NeutrinoTrackingEff_module.cc.

TH1D* DUNE::NeutrinoTrackingEff::h_kaon_length
private

Definition at line 161 of file NeutrinoTrackingEff_module.cc.

TH1D* DUNE::NeutrinoTrackingEff::h_kaonwtrk_length
private

Definition at line 163 of file NeutrinoTrackingEff_module.cc.

TH1D* DUNE::NeutrinoTrackingEff::h_michel_length
private

Definition at line 162 of file NeutrinoTrackingEff_module.cc.

TH1D* DUNE::NeutrinoTrackingEff::h_michelwtrk_length
private

Definition at line 164 of file NeutrinoTrackingEff_module.cc.

TH1D* DUNE::NeutrinoTrackingEff::h_muon_length
private

Definition at line 129 of file NeutrinoTrackingEff_module.cc.

TH1D* DUNE::NeutrinoTrackingEff::h_muonwtrk_length
private

Definition at line 133 of file NeutrinoTrackingEff_module.cc.

TH1D* DUNE::NeutrinoTrackingEff::h_pionm_length
private

Definition at line 132 of file NeutrinoTrackingEff_module.cc.

TH1D* DUNE::NeutrinoTrackingEff::h_pionmwtrk_length
private

Definition at line 136 of file NeutrinoTrackingEff_module.cc.

TH1D* DUNE::NeutrinoTrackingEff::h_pionp_length
private

Definition at line 131 of file NeutrinoTrackingEff_module.cc.

TH1D* DUNE::NeutrinoTrackingEff::h_pionpwtrk_length
private

Definition at line 135 of file NeutrinoTrackingEff_module.cc.

TH1D* DUNE::NeutrinoTrackingEff::h_Pkaon_den
private

Definition at line 151 of file NeutrinoTrackingEff_module.cc.

TH1D* DUNE::NeutrinoTrackingEff::h_Pkaon_num
private

Definition at line 152 of file NeutrinoTrackingEff_module.cc.

TH1D* DUNE::NeutrinoTrackingEff::h_Pmichel_e_den
private

Definition at line 153 of file NeutrinoTrackingEff_module.cc.

TH1D* DUNE::NeutrinoTrackingEff::h_Pmichel_e_num
private

Definition at line 154 of file NeutrinoTrackingEff_module.cc.

TH1D* DUNE::NeutrinoTrackingEff::h_Pmu_den
private

Definition at line 105 of file NeutrinoTrackingEff_module.cc.

TH1D* DUNE::NeutrinoTrackingEff::h_Pmu_num
private

Definition at line 106 of file NeutrinoTrackingEff_module.cc.

TH1D* DUNE::NeutrinoTrackingEff::h_Ppion_minus_den
private

Definition at line 113 of file NeutrinoTrackingEff_module.cc.

TH1D* DUNE::NeutrinoTrackingEff::h_Ppion_minus_num
private

Definition at line 114 of file NeutrinoTrackingEff_module.cc.

TH1D* DUNE::NeutrinoTrackingEff::h_Ppion_plus_den
private

Definition at line 111 of file NeutrinoTrackingEff_module.cc.

TH1D* DUNE::NeutrinoTrackingEff::h_Ppion_plus_num
private

Definition at line 112 of file NeutrinoTrackingEff_module.cc.

TH1D* DUNE::NeutrinoTrackingEff::h_Pproton_den
private

Definition at line 109 of file NeutrinoTrackingEff_module.cc.

TH1D* DUNE::NeutrinoTrackingEff::h_Pproton_num
private

Definition at line 110 of file NeutrinoTrackingEff_module.cc.

TH1D* DUNE::NeutrinoTrackingEff::h_proton_length
private

Definition at line 130 of file NeutrinoTrackingEff_module.cc.

TH1D* DUNE::NeutrinoTrackingEff::h_protonwtrk_length
private

Definition at line 134 of file NeutrinoTrackingEff_module.cc.

TH1D* DUNE::NeutrinoTrackingEff::h_theta_den
private

Definition at line 107 of file NeutrinoTrackingEff_module.cc.

TH1D* DUNE::NeutrinoTrackingEff::h_theta_num
private

Definition at line 108 of file NeutrinoTrackingEff_module.cc.

TH1D* DUNE::NeutrinoTrackingEff::h_trackRes_kaon
private

Definition at line 157 of file NeutrinoTrackingEff_module.cc.

TH1D* DUNE::NeutrinoTrackingEff::h_trackRes_lepton
private

Definition at line 124 of file NeutrinoTrackingEff_module.cc.

TH1D* DUNE::NeutrinoTrackingEff::h_trackRes_michel
private

Definition at line 160 of file NeutrinoTrackingEff_module.cc.

TH1D* DUNE::NeutrinoTrackingEff::h_trackRes_pion_minus
private

Definition at line 127 of file NeutrinoTrackingEff_module.cc.

TH1D* DUNE::NeutrinoTrackingEff::h_trackRes_pion_plus
private

Definition at line 126 of file NeutrinoTrackingEff_module.cc.

TH1D* DUNE::NeutrinoTrackingEff::h_trackRes_proton
private

Definition at line 125 of file NeutrinoTrackingEff_module.cc.

double DUNE::NeutrinoTrackingEff::MC_incoming_P[4]
private

Definition at line 85 of file NeutrinoTrackingEff_module.cc.

int DUNE::NeutrinoTrackingEff::MC_incoming_PDG
private

Definition at line 84 of file NeutrinoTrackingEff_module.cc.

int DUNE::NeutrinoTrackingEff::MC_isCC
private

Definition at line 83 of file NeutrinoTrackingEff_module.cc.

int DUNE::NeutrinoTrackingEff::MC_kaonID
private

Definition at line 93 of file NeutrinoTrackingEff_module.cc.

double DUNE::NeutrinoTrackingEff::MC_kaonP
private

Definition at line 100 of file NeutrinoTrackingEff_module.cc.

int DUNE::NeutrinoTrackingEff::MC_leading_PionMinusID
private

Definition at line 91 of file NeutrinoTrackingEff_module.cc.

double DUNE::NeutrinoTrackingEff::MC_leading_PionMinusP
private

Definition at line 99 of file NeutrinoTrackingEff_module.cc.

int DUNE::NeutrinoTrackingEff::MC_leading_PionPlusID
private

Definition at line 90 of file NeutrinoTrackingEff_module.cc.

double DUNE::NeutrinoTrackingEff::MC_leading_PionPlusP
private

Definition at line 97 of file NeutrinoTrackingEff_module.cc.

int DUNE::NeutrinoTrackingEff::MC_leading_protonID
private

Definition at line 89 of file NeutrinoTrackingEff_module.cc.

double DUNE::NeutrinoTrackingEff::MC_leading_ProtonP
private

Definition at line 98 of file NeutrinoTrackingEff_module.cc.

double DUNE::NeutrinoTrackingEff::MC_lepton_startMomentum[4]
private

Definition at line 87 of file NeutrinoTrackingEff_module.cc.

int DUNE::NeutrinoTrackingEff::MC_leptonID
private

Definition at line 92 of file NeutrinoTrackingEff_module.cc.

double DUNE::NeutrinoTrackingEff::MC_leptonP
private

Definition at line 96 of file NeutrinoTrackingEff_module.cc.

int DUNE::NeutrinoTrackingEff::MC_michelID
private

Definition at line 94 of file NeutrinoTrackingEff_module.cc.

double DUNE::NeutrinoTrackingEff::MC_michelP
private

Definition at line 101 of file NeutrinoTrackingEff_module.cc.

double DUNE::NeutrinoTrackingEff::MC_vertex[4]
private

Definition at line 86 of file NeutrinoTrackingEff_module.cc.


The documentation for this class was generated from the following file: