MuonTrackingEff_module.cc
Go to the documentation of this file.
1 //
2 // **Muon Tracking Efficiency module**
3 // This module is based on NeutrinoTrackingEff_module.cc from A. Higuera.
4 // I have changed the module so that it only calculates the Muon Tracking
5 // Efficiency by checking the Completeness and Purity (see function: "truthMatcher")
6 // and track length (see function: "truthLength") of the leading reconstructed muon
7 // track (This is the one with the highest Completeness).
8 // In case the leading muon track failed and there is more than one reconstructed
9 // muon track (e.g. due tue to track splitting in the reconstruction bc of a kink
10 // after multiple scattering), I check the efficiency criteria for the sum of the
11 // leading and second reconstructed muon track.
12 //
13 // Christoph Alt
14 // christoph.alt@cern.ch
15 
16 // LArSoft includes
23 
24 // Framework includes
30 #include "art_root_io/TFileService.h"
31 #include "canvas/Persistency/Common/FindManyP.h"
32 #include "fhiclcpp/ParameterSet.h"
34 
35 // ROOT includes
36 #include "TH1.h"
37 #include "TH2.h"
38 #include "TStyle.h"
39 #include "TVector3.h"
40 
41 #define MAX_TRACKS 1000
42 using namespace std;
43 
44 //========================================================================
45 
46 namespace DUNE {
47 
49  public:
50  explicit MuonTrackingEff(fhicl::ParameterSet const& pset);
51 
52  private:
53  void beginJob() override;
54  void endJob() override;
55  void beginRun(const art::Run& run) override;
56  void analyze(const art::Event& evt) override;
57 
58  void processEff(const art::Event& evt, bool& isFiducial);
59 
60  void truthMatcher(detinfo::DetectorClocksData const& clockData,
63  const simb::MCParticle*& MCparticle,
64  double& Purity,
65  double& Completeness,
66  double& TotalRecoEnergy);
67 
68  void FuncDistanceAndAngleBetweenTracks(art::Ptr<recob::Track> Track1,
70  double& TempDistanceBetweenTracks,
71  double& TempAngleBetweenTracks,
72  double& TempCriteriaTwoTracks);
73 
74  void FuncDistanceAndAngleBetweenTruthAndRecoTrack(const simb::MCParticle*& MCparticle,
76  double& TempDistanceBetweenTruthAndRecoTrack,
77  double& TempAngleBeetweenTruthAndRecoTrack);
78 
79  double truthLength(const simb::MCParticle* MCparticle);
80 
81  bool insideFV(double vertex[4]);
82 
83  void doEfficiencies();
84 
85  // the parameters we'll read from the .fcl
89 
90  // int MC_isCC;
91  // int MC_incoming_PDG;
92  // double MC_incoming_P[4];
93  double MCTruthMuonVertex[4];
94  // double MCTruthMuonStartMomentum[4];
95 
98 
99  // double MCTruthMuonMomentumXZ=0;
100  // double MCTruthMuonMomentumYZ=0;
101  double MCTruthMuonThetaXZ = 0;
102  double MCTruthMuonThetaYZ = 0;
103 
104  //Counts to calculate efficiency and failed criteria
105  int EventCounter = 0;
106 
107  int CountMCTruthMuon = 0;
108  int CountRecoMuon = 0;
109 
110  int CountGoodLeadingMuonTrack = 0;
111  int CountNoRecoTracks = 0;
112  int CountNoMuonTracks = 0;
113  int CountBadLeadingMuonTrack = 0;
114  int CountCompleteness = 0;
115  int CountPurity = 0;
116  int CountTrackLengthTooShort = 0;
117  int CountTrackLengthTooLong = 0;
118 
119  int CountBadLeadingMuonTrackAndOnlyOneMuonTrack = 0;
120  int CountBadLeadingMuonTrackButLeadingPlusSecondGood = 0;
121  int CountBadLeadingMuonTrackAndLeadingPlusSecondBad = 0;
122 
123  int CountBadLeadingMuonTrackAndLeadingPlusSecondBadCompleteness = 0;
124  int CountBadLeadingMuonTrackAndLeadingPlusSecondBadPurity = 0;
125  int CountBadLeadingMuonTrackAndLeadingPlusSecondBadTrackTooShort = 0;
126  int CountBadLeadingMuonTrackAndLeadingPlusSecondBadTrackTooLong = 0;
127 
128  int CountBadLeadingMuonTrackAndOnlyOneMuonTrackCompleteness = 0;
129  int CountBadLeadingMuonTrackAndOnlyOneMuonTrackPurity = 0;
130  int CountBadLeadingMuonTrackAndOnlyOneMuonTrackTrackTooShort = 0;
131  int CountBadLeadingMuonTrackAndOnlyOneMuonTrackTrackTooLong = 0;
132 
133  double Criteria;
134  // int NMuonTracksTooShort=0;
135 
136  int GoodEvents1MuonTrack = 0;
137  int GoodEvents2MuonTrack = 0;
138  int GoodEvents3MuonTrack = 0;
139  int GoodEvents4OrMoreMuonTrack = 0;
140 
141  int BadEvents0MuonTrack = 0;
142  int BadEvents1MuonTrack = 0;
143  int BadEvents2MuonTrack = 0;
144  int BadEvents3MuonTrack = 0;
145  int BadEvents4OrMoreMuonTrack = 0;
146 
147  //TH1Ds
148  //Single Criteria and total reco energy
149  TH1D* h_Purity;
151  TH1D* h_TrackRes;
154  TH1D* h_VertexRes;
156 
157  //Efficiency ThetaXZ
161 
162  //Efficiency ThetaYZ
166 
167  //Efficiency SinThetaYZ
171 
172  //TH2Ds
173  //Efficiency ThetaXZ vs. ThetaYZ
178 
179  //Efficiency ThetaXZ vs. SinThetaYZ
184 
185  //Efficiency ThetaXZ vs. ThetaYZ after summing up leading plus second
189 
190  //Difference in efficiency before and after summing up leading plus second: ThetaXZ vs. ThetaYZ
192 
193  //Efficiency ThetaXZ vs. SinThetaYZ after summing up leading plus second
197 
198  //Failed Criteria
205 
206  //Criteria vs. NRecoTrack
210 
211  //Criteria vs. NMuonTrack
215 
216  //NoMuonTrack: Max length of no muon track vs. PDG code of that track (MC truth)
218 
219  //Stitching variables: all events
226 
227  //Stitching variables: bad events
235 
236  //Stitching variables: good events
241 
242  //Stitching variables: bad events but leading plus second ok
244  TH2D*
247 
248  //Stitching variables: bad events and leading plus second not ok
252 
253  float fFidVolCutX;
254  float fFidVolCutY;
255  float fFidVolCutZ;
256 
257  float fFidVolXmin;
258  float fFidVolXmax;
259  float fFidVolYmin;
260  float fFidVolYmax;
261  float fFidVolZmin;
262  float fFidVolZmax;
263 
265 
266  //My histograms
267  int NThetaXZBins = 36;
268  int ThetaXZBinMin = 0;
269  int ThetaXZBinMax = 360;
270 
271  int NThetaYZBins = 18;
272  int ThetaYZBinMin = -90;
273  int ThetaYZBinMax = 90;
274 
275  int NSinThetaYZBins = 18;
276  int SinThetaYZBinMin = -1;
277  int SinThetaYZBinMax = 1;
278 
279  int NCriteriaBins = 13;
280  double CriteriaBinMin = -0.25;
281  double CriteriaBinMax = 6.25;
282 
283  int NRecoTracksBins = 19;
284  double RecoTracksBinMin = -0.25;
285  double RecoTracksBinMax = 9.25;
286 
287  }; // class MuonTrackingEff
288 
289  //========================================================================
290  MuonTrackingEff::MuonTrackingEff(fhicl::ParameterSet const& p) : EDAnalyzer(p)
291  {
292  fMCTruthModuleLabel = p.get<std::string>("MCTruthModuleLabel");
293  fTrackModuleLabel = p.get<std::string>("TrackModuleLabel");
294  fMuonPDGCode = p.get<int>("MuonPDGCode");
295  fFidVolCutX = p.get<float>("FidVolCutX");
296  fFidVolCutY = p.get<float>("FidVolCutY");
297  fFidVolCutZ = p.get<float>("FidVolCutZ");
298  }
299  //========================================================================
300  void
302  {
303  std::cout << "job begin..." << std::endl;
304  // Get geometry.
305  auto const* geo = lar::providerFrom<geo::Geometry>();
306  // Define histogram boundaries (cm).
307  // For now only draw cryostat=0.
308  double minx = 1e9;
309  double maxx = -1e9;
310  double miny = 1e9;
311  double maxy = -1e9;
312  double minz = 1e9;
313  double maxz = -1e9;
314  for (size_t i = 0; i < geo->NTPC(); ++i) {
315  double local[3] = {0., 0., 0.};
316  double world[3] = {0., 0., 0.};
317  const geo::TPCGeo& tpc = geo->TPC(i);
318  tpc.LocalToWorld(local, world);
319  if (minx > world[0] - geo->DetHalfWidth(i)) minx = world[0] - geo->DetHalfWidth(i);
320  if (maxx < world[0] + geo->DetHalfWidth(i)) maxx = world[0] + geo->DetHalfWidth(i);
321  if (miny > world[1] - geo->DetHalfHeight(i)) miny = world[1] - geo->DetHalfHeight(i);
322  if (maxy < world[1] + geo->DetHalfHeight(i)) maxy = world[1] + geo->DetHalfHeight(i);
323  if (minz > world[2] - geo->DetLength(i) / 2.) minz = world[2] - geo->DetLength(i) / 2.;
324  if (maxz < world[2] + geo->DetLength(i) / 2.) maxz = world[2] + geo->DetLength(i) / 2.;
325  }
326 
327  fFidVolXmin = minx + fFidVolCutX;
328  fFidVolXmax = maxx - fFidVolCutX;
329  fFidVolYmin = miny + fFidVolCutY;
330  fFidVolYmax = maxy - fFidVolCutY;
331  fFidVolZmin = minz + fFidVolCutZ;
332  fFidVolZmax = maxz - fFidVolCutZ;
333 
334  std::cout << "Fiducial volume:"
335  << "\n"
336  << fFidVolXmin << "\t< x <\t" << fFidVolXmax << "\n"
337  << fFidVolYmin << "\t< y <\t" << fFidVolYmax << "\n"
338  << fFidVolZmin << "\t< z <\t" << fFidVolZmax << "\n";
339 
341 
342  //TH1D's
343  gStyle->SetTitleOffset(1.3, "Y");
344 
345  //Single Criteria and total reco energy
346  h_Purity =
347  tfs->make<TH1D>("h_Purity", "All events: Purity vs. # events; Purity; # events", 60, 0, 1.2);
348 
349  h_Completeness = tfs->make<TH1D>(
350  "h_Completeness", "All events: Completeness vs # events; Completeness; # events", 60, 0, 1.2);
351  h_Completeness->SetLineColor(kBlue);
352 
353  h_TrackRes =
354  tfs->make<TH1D>("h_TrackRes",
355  "All events: L_{reco}/L_{truth} vs. # events; L_{reco}/L_{truth}; # events;",
356  75,
357  0,
358  1.5);
359  h_TrackRes->SetLineColor(kRed);
360 
361  h_TotalRecoEnergy = tfs->make<TH1D>("h_TotalRecoEnergy",
362  "All events: Total reco energy (sum of all hits in all "
363  "tracks) vs. # events; E_{reco., tot.} [MeV]; # events",
364  100,
365  0,
366  1000);
367 
368  h_TruthLength = tfs->make<TH1D>(
369  "h_TruthLength",
370  "All events: truth muon length vs. # events; truth muon length [cm]; # events",
371  100,
372  0,
373  300);
374 
375  h_VertexRes = tfs->make<TH1D>(
376  "h_VertexRes",
377  "All events: Vertex residuals vs. # events; #Delta vertex_{truth-teco} [cm]; # events",
378  300,
379  0,
380  300);
381 
382  h_DirectionRes = tfs->make<TH1D>(
383  "h_DirectionRes",
384  "All events: Angular residuals vs. # events; #Delta#theta_{truth-reco} [#circ]; # events",
385  180,
386  0,
387  180);
388 
389  //Efficiency ThetaXZ
391  tfs->make<TH1D>("h_Efficiency_ThetaXZ",
392  "Muon reco efficiency vs. #theta_{XZ}; #theta_{XZ} [#circ]; Efficiency",
393  NThetaXZBins,
395  ThetaXZBinMax);
396  h_ThetaXZ_den =
397  tfs->make<TH1D>("h_ThetaXZ_den",
398  "# generated muons vs. #theta_{XZ}; #theta_{XZ} [#circ]; # generated muons",
399  NThetaXZBins,
401  ThetaXZBinMax);
402  h_ThetaXZ_num = tfs->make<TH1D>(
403  "h_ThetaXZ_num",
404  "# reconstructed muons vs. #theta_{XZ}; #theta_{XZ} [#circ]; # reconstructed muons",
405  NThetaXZBins,
407  ThetaXZBinMax);
408 
409  //Efficiency ThetaYZ
410  h_Efficiency_ThetaYZ = tfs->make<TH1D>(
411  "h_Efficiency_ThetaYZ",
412  "Muon reco efficiency vs. #theta_{YZ}; #theta_{YZ} [#circ]; Muon reco efficiency",
413  NThetaYZBins,
415  ThetaYZBinMax);
416  ;
417  h_ThetaYZ_den =
418  tfs->make<TH1D>("h_ThetaYZ_den",
419  "# generated muons vs. #theta_{YZ}; #theta_{YZ} [#circ]; # generated muons",
420  NThetaYZBins,
422  ThetaYZBinMax);
423  h_ThetaYZ_num = tfs->make<TH1D>(
424  "h_ThetaYZ_num",
425  "# reconstructed muons vs. #theta_{YZ}; #theta_{YZ} [#circ]; # reconstructed muons",
426  NThetaYZBins,
428  ThetaYZBinMax);
429 
430  //Efficiency SinThetaYZ
431  h_Efficiency_SinThetaYZ = tfs->make<TH1D>(
432  "h_Efficiency_SinThetaYZ",
433  "Muon reco efficiency vs. sin(#theta_{YZ}); sin(#theta_{YZ}); Muon reco efficiency",
437  ;
439  tfs->make<TH1D>("h_SinThetaYZ_den",
440  "# generated muons vs. sin(#theta_{YZ}); sin(#theta_{YZ}); # generated muons",
444  h_SinThetaYZ_num = tfs->make<TH1D>(
445  "h_SinThetaYZ_num",
446  "# reconstructed muons vs. sin(#theta_{YZ}); sin(#theta_{YZ}); # reconstructed muons",
450 
451  h_Purity->Sumw2();
452  h_Completeness->Sumw2();
453  h_TrackRes->Sumw2();
454  h_TotalRecoEnergy->Sumw2();
455  h_TruthLength->Sumw2();
456  h_VertexRes->Sumw2();
457  h_DirectionRes->Sumw2();
458 
459  h_Efficiency_SinThetaYZ->Sumw2();
460  h_SinThetaYZ_den->Sumw2();
461  h_SinThetaYZ_num->Sumw2();
462 
463  h_Efficiency_ThetaXZ->Sumw2();
464  h_ThetaXZ_den->Sumw2();
465  h_ThetaXZ_num->Sumw2();
466 
467  h_Efficiency_ThetaYZ->Sumw2();
468  h_ThetaYZ_den->Sumw2();
469  h_ThetaYZ_num->Sumw2();
470 
471  //TH2D's
472  //Efficiency and Failed Reconstruction ThetaXZ vs. ThetaYZ
473  h_Efficiency_ThetaXZ_ThetaYZ = tfs->make<TH2D>(
474  "h_Efficiency_ThetaXZ_ThetaYZ",
475  "Muon reco efficiency: #theta_{XZ} vs. #theta_{YZ}; #theta_{XZ} [#circ]; #theta_{YZ} [#circ]",
476  NThetaXZBins,
479  NThetaYZBins,
481  ThetaYZBinMax);
482  h_Efficiency_ThetaXZ_ThetaYZ->SetOption("colz");
483 
484  h_ThetaXZ_ThetaYZ_den = tfs->make<TH2D>(
485  "h_ThetaXZ_ThetaYZ_den",
486  "# generated muons: #theta_{XZ} vs. #theta_{YZ}; #theta_{XZ} [#circ]; #theta_{YZ} [#circ]",
487  NThetaXZBins,
490  NThetaYZBins,
492  ThetaYZBinMax);
493  h_ThetaXZ_ThetaYZ_den->SetOption("colz");
494 
495  h_ThetaXZ_ThetaYZ_num = tfs->make<TH2D>("h_ThetaXZ_ThetaYZ_num",
496  "# reconstructed muons: #theta_{XZ} vs. #theta_{YZ}; "
497  "#theta_{XZ} [#circ]; #theta_{YZ} [#circ]",
498  NThetaXZBins,
501  NThetaYZBins,
503  ThetaYZBinMax);
504  h_ThetaXZ_ThetaYZ_num->SetOption("colz");
505 
507  tfs->make<TH2D>("h_FailedReconstruction_ThetaXZ_ThetaYZ",
508  "# failed reconstructions: #theta_{XZ} vs. #theta_{YZ}; #theta_{XZ} [#circ]; "
509  "#theta_{YZ} [#circ]",
510  NThetaXZBins,
513  NThetaYZBins,
515  ThetaYZBinMax);
516  h_FailedReconstruction_ThetaXZ_ThetaYZ->SetOption("colz");
517 
518  //Efficiency and Failed Reconstruction ThetaXZ vs. SinThetaYZ
520  tfs->make<TH2D>("h_Efficiency_ThetaXZ_SinThetaYZ",
521  "Muon reco efficiency: #theta_{XZ} vs. sin(#theta_{YZ}); #theta_{XZ} "
522  "[#circ]; sin(#theta_{YZ})",
523  NThetaXZBins,
529  h_Efficiency_ThetaXZ_SinThetaYZ->SetOption("colz");
530 
531  h_ThetaXZ_SinThetaYZ_den = tfs->make<TH2D>(
532  "h_ThetaXZ_SinThetaYZ_den",
533  "# generated muons: #theta_{XZ} vs. sin(#theta_{YZ}); #theta_{XZ} [#circ]; sin(#theta_{YZ})",
534  NThetaXZBins,
540  h_ThetaXZ_SinThetaYZ_den->SetOption("colz");
541 
543  tfs->make<TH2D>("h_ThetaXZ_SinThetaYZ_num",
544  "# reconstructed muons: #theta_{XZ} vs. sin(#theta_{YZ}); #theta_{XZ} "
545  "[#circ]; sin(#theta_{YZ}); #theta_{XZ} [#circ]; sin(#theta_{YZ})",
546  NThetaXZBins,
552  h_ThetaXZ_SinThetaYZ_num->SetOption("colz");
553 
555  tfs->make<TH2D>("h_FailedReconstruction_ThetaXZ_SinThetaYZ",
556  "# failed reconstructions: #theta_{XZ} vs. sin(#theta_{YZ}); #theta_{XZ} "
557  "[#circ]; sin(#theta_{YZ})",
558  NThetaXZBins,
565 
566  //Efficiency ThetaXZ vs. ThetaYZ after summing up leading plus second
568  tfs->make<TH2D>("h_Efficiency_ThetaXZ_ThetaYZ_LeadingPlusSecond",
569  "Muon reco efficiency after stitching: #theta_{XZ} vs. #theta_{YZ}; "
570  "#theta_{XZ} [#circ]; #theta_{YZ} [#circ]",
571  NThetaXZBins,
574  NThetaYZBins,
576  ThetaYZBinMax);
578 
580  tfs->make<TH2D>("h_ThetaXZ_ThetaYZ_LeadingPlusSecondOk",
581  "# reconstructed muons after stitching (failed before stitching): "
582  "#theta_{XZ} vs #theta_{YZ}; #theta_{XZ} [#circ]; #theta_{YZ} [#circ]",
583  NThetaXZBins,
586  NThetaYZBins,
588  ThetaYZBinMax);
589  h_ThetaXZ_ThetaYZ_LeadingPlusSecondOk->SetOption("colz");
590 
592  tfs->make<TH2D>("h_ThetaXZ_ThetaYZ_LeadingPlusSecondOk_num",
593  "# reconstructed muons after stitching: #theta_{XZ} vs. #theta_{YZ}; "
594  "#theta_{XZ} [#circ]; #theta_{YZ} [#circ]",
595  NThetaXZBins,
598  NThetaYZBins,
600  ThetaYZBinMax);
602 
603  //Efficiency ThetaXZ vs. SinThetaYZ after summing up leading plus second
605  tfs->make<TH2D>("h_Efficiency_ThetaXZ_SinThetaYZ_LeadingPlusSecond",
606  "Muon reco efficiency after stitching: #theta_{XZ} vs. sin(#theta_{YZ}); "
607  "#theta_{XZ} [#circ]; sin(#theta_{YZ})",
608  NThetaXZBins,
615 
617  tfs->make<TH2D>("h_ThetaXZ_SinThetaYZ_LeadingPlusSecondOk",
618  "# reconstructed muons after stitching (failed before stitching): "
619  "#theta_{XZ} vs. sin(#theta_{YZ}); #theta_{XZ} [#circ]; sin(#theta_{YZ})",
620  NThetaXZBins,
626  h_ThetaXZ_SinThetaYZ_LeadingPlusSecondOk->SetOption("colz");
627 
629  tfs->make<TH2D>("h_ThetaXZ_SinThetaYZ_LeadingPlusSecondOk_num",
630  "# reconstructed muons after stitching: #theta_{XZ} vs. sin(#theta_{YZ}); "
631  "#theta_{XZ} [#circ]; sin(#theta_{YZ})",
632  NThetaXZBins,
639 
640  //Difference in efficiency before and after summing up leading plus second: ThetaXZ vs. ThetaYZ
642  tfs->make<TH2D>("h_Efficiency_ThetaXZ_ThetaYZ_DifferenceLeadingAndLeadingPlusSecond",
643  "Muon reco efficiency: difference before and after stitching: #theta_{XZ} "
644  "vs. #theta_{YZ}; #theta_{XZ} [#circ]; #theta_{YZ} [#circ]",
645  NThetaXZBins,
648  NThetaYZBins,
650  ThetaYZBinMax);
652 
653  //Failed Criteria
655  tfs->make<TH2D>("h_NoRecoTrackAtAll_ThetaXZ_SinThetaYZ",
656  "# events with no reco track at all: #theta_{XZ} vs. sin(#theta_{YZ}); "
657  "#theta_{XZ} [#circ]; sin(#theta_{YZ})",
658  NThetaXZBins,
664  h_NoRecoTrackAtAll_ThetaXZ_SinThetaYZ->SetOption("colz");
665 
667  tfs->make<TH2D>("h_NoMuonTrack_ThetaXZ_SinThetaYZ",
668  "# events with no muon track: #theta_{XZ} vs. sin(#theta_{YZ}); #theta_{XZ} "
669  "[#circ]; sin(#theta_{YZ})",
670  NThetaXZBins,
676  h_NoMuonTrack_ThetaXZ_SinThetaYZ->SetOption("colz");
677 
679  tfs->make<TH2D>("h_TrackTooShort_ThetaXZ_SinThetaYZ",
680  "# events with L_{reco}/L_{truth} < 75%: #theta_{XZ} vs. sin(#theta_{YZ}); "
681  "#theta_{XZ} [#circ]; sin(#theta_{YZ})",
682  NThetaXZBins,
688  h_TrackTooShort_ThetaXZ_SinThetaYZ->SetOption("colz");
689 
691  tfs->make<TH2D>("h_TrackTooLong_ThetaXZ_SinThetaYZ",
692  "#events with L_{reco}/L_{truth} > 125%: #theta_{XZ} vs. sin(#theta_{YZ}); "
693  "#theta_{XZ} [#circ]; sin(#theta_{YZ})",
694  NThetaXZBins,
700  h_TrackTooLong_ThetaXZ_SinThetaYZ->SetOption("colz");
701 
703  tfs->make<TH2D>("h_Completeness_ThetaXZ_SinThetaYZ",
704  "# events with Completeness < 50%: #theta_{XZ} vs. sin(#theta_{YZ}); "
705  "#theta_{XZ} [#circ]; sin(#theta_{YZ})",
706  NThetaXZBins,
712  h_Completeness_ThetaXZ_SinThetaYZ->SetOption("colz");
713 
715  tfs->make<TH2D>("h_Purity_ThetaXZ_SinThetaYZ",
716  "# events with Purity < 50%: #theta_{XZ} vs. sin(#theta_{YZ}); #theta_{XZ} "
717  "[#circ]; sin(#theta_{YZ})",
718  NThetaXZBins,
724  h_Purity_ThetaXZ_SinThetaYZ->SetOption("colz");
725 
726  //Criteria vs. NRecoTrack
728  tfs->make<TH2D>("h_Criteria_NRecoTrack",
729  "Ratio: criteria vs. # reco tracks; Criteria; # reco tracks",
736  h_Criteria_NRecoTrack->SetOption("colz");
737 
739  tfs->make<TH2D>("h_Criteria_NRecoTrack_num",
740  "# events: criteria vs. # reco tracks; Criteria; # reco tracks",
747  h_Criteria_NRecoTrack_num->SetOption("colz");
748 
749  h_Criteria_NRecoTrack_den = tfs->make<TH2D>("h_Criteria_NRecoTrack_den",
750  "Divider histogram; Criteria; # reco tracks",
757  h_Criteria_NRecoTrack_den->SetOption("colz");
758 
759  //Criteria vs. NMuonTrack
761  tfs->make<TH2D>("h_Criteria_NMuonTrack",
762  "Ratio: criteria vs. # muon tracks; Criteria; # muon tracks",
769  h_Criteria_NMuonTrack->SetOption("colz");
770 
772  tfs->make<TH2D>("h_Criteria_NMuonTrack_num",
773  "# events: criteria vs. # muon tracks; Criteria; # muon tracks",
780  h_Criteria_NMuonTrack_num->SetOption("colz");
781 
782  h_Criteria_NMuonTrack_den = tfs->make<TH2D>("h_Criteria_NMuonTrack_den",
783  "Divider histogram; Criteria; # muon tracks",
790  h_Criteria_NMuonTrack_den->SetOption("colz");
791 
792  //NoMuonTrack: Max length of no muon track vs. PDG code of that track (MC truth)
793  h_NoMuonTrack_MaxTrackLength_PDGCode = tfs->make<TH2D>(
794  "h_NoMuonTrack_MaxTrackLength_PDGCode",
795  "Events with no muon track: L_{reco, max} vs. PDG Code; L_{reco} [cm]; PDG Code",
796  100,
797  0.,
798  200.,
799  200,
800  -100.,
801  100.);
802  h_NoMuonTrack_MaxTrackLength_PDGCode->SetOption("colz");
803 
804  //Stitching variables: all events
806  tfs->make<TH2D>("h_MuonTrackStitching_TrackRes_Completeness",
807  "All events: L_{reco}/L_{truth} (leading) vs. Completeness (leading); "
808  "L_{reco}/L_{truth} (leading); Completeness (leading)",
809  150,
810  0.,
811  1.5,
812  150,
813  0.,
814  1.5);
816 
818  tfs->make<TH2D>("h_MuonTrackStitching_TrackResLeading_TrackResSecond",
819  "All events: L_{reco}/L_{truth} (leading) vs. L_{reco}/L_{truth} (second); "
820  "L_{reco}/L_{truth} (leading); L_{reco}/L_{truth} (second)",
821  150,
822  0.,
823  1.5,
824  150,
825  0.,
826  1.5);
828 
830  tfs->make<TH2D>("h_MuonTrackStitching_Distance_Angle",
831  "All events: distance vs. angle b/w leading and second muon track; Distance "
832  "[cm]; Angle [#circ]",
833  100,
834  0.,
835  100.,
836  100,
837  0.,
838  180.);
839  h_MuonTrackStitching_Distance_Angle->SetOption("colz");
840 
842  tfs->make<TH2D>("h_MuonTrackStitching_TrackResSecondMuon_Angle",
843  "All events: L_{reco}/L_{truth} (second) vs. angle; L_{reco}/L_{truth} "
844  "(second); Angle [#circ]",
845  150,
846  0.,
847  1.5,
848  180,
849  0,
850  180.);
852 
854  "h_MuonTrackStitching_CompletenessSecondMuon_Angle",
855  "All events: Completeness (second) vs. angle; Completeness (second); Angle [#circ]",
856  120,
857  0.,
858  1.2,
859  180,
860  0,
861  180.);
863 
865  tfs->make<TH2D>("h_MuonTrackStitching_CriteriaTwoTracks_Angle",
866  "All events: CriteriaTwoTracks vs. angle; Criteria; Angle [#circ]",
867  7,
868  0.75,
869  4.25,
870  180,
871  0,
872  180.);
874 
875  //Stitching variables: bad events
877  tfs->make<TH2D>("h_MuonTrackStitching_FailedCriteria_TrackRes_Completeness",
878  "Bad events: L_{reco}/L_{truth} (leading) vs. Completeness (leading); "
879  "L_{reco}/L_{truth} (leading); Completeness (leading)",
880  150,
881  0.,
882  1.5,
883  150,
884  0.,
885  1.5);
887 
889  tfs->make<TH2D>("h_MuonTrackStitching_FailedCriteria_Distance_Angle",
890  "Bad events: distance vs. angle b/w leading and second muon track; Distance "
891  "[cm]; Angle [#circ]",
892  100,
893  0.,
894  100.,
895  100,
896  0.,
897  180.);
899 
901  tfs->make<TH2D>("h_MuonTrackStitching_FailedCriteria_TrackResLeading_TrackResSecond",
902  "Bad events: L_{reco}/L_{truth} (leading) vs. L_{reco}/L_{truth} (second); "
903  "L_{reco}/L_{truth} (leading); L_{reco}/L_{truth} (second)",
904  150,
905  0.,
906  1.5,
907  150,
908  0.,
909  1.5);
911 
913  tfs->make<TH2D>("*h_MuonTrackStitching_FailedCriteria_CompletenessLeading_CompletenessSecond",
914  "Bad events: Completeness (leading) vs. Completeness (second); Completeness "
915  "(leading); Completeness (second)",
916  150,
917  0.,
918  1.5,
919  150,
920  0.,
921  1.5);
923 
925  tfs->make<TH2D>("h_MuonTrackStitching_FailedCriteria_TrackResSecondMuon_Angle",
926  "Bad events: L_{reco}/L_{truth} (second) vs. angle; L_{reco}/L_{truth} "
927  "(second); Angle [#circ]",
928  150,
929  0.,
930  1.5,
931  180,
932  0,
933  180.);
935 
937  "h_MuonTrackStitching_FailedCriteria_CompletenessSecondMuon_Angle",
938  "Bad events: Completeness (second) vs. angle; Completeness (second); Angle [#circ]",
939  120,
940  0.,
941  1.2,
942  180,
943  0,
944  180.);
946 
948  tfs->make<TH2D>("h_MuonTrackStitching_FailedCriteria_CriteriaTwoTracks_Angle",
949  "Bad events: CriteriaTwoTracks vs. angle; Criteria; Angle [#circ]",
950  7,
951  0.75,
952  4.25,
953  180,
954  0,
955  180.);
957 
958  //Stitching variables: good events
960  tfs->make<TH2D>("h_MuonTrackStitching_MatchedCriteria_TrackRes_Completeness",
961  "Good events: L_{reco}/L_{truth} (leading) vs. Completeness (leading); "
962  "L_{reco}/L_{truth} (leading); Completeness (leading)",
963  150,
964  0.,
965  1.5,
966  150,
967  0.,
968  1.5);
970 
972  tfs->make<TH2D>("h_MuonTrackStitching_MatchedCriteria_Distance_Angle",
973  "Good events: distance vs. angle b/w leading and second muon track; Distance "
974  "[cm]; Angle [#circ]",
975  100,
976  0.,
977  100.,
978  100,
979  0.,
980  180.);
982 
984  tfs->make<TH2D>("h_MuonTrackStitching_MatchedCriteria_TrackResLeading_TrackResSecond",
985  "Good events: L_{reco}/L_{truth} (leading) vs. L_{reco}/L_{truth} (second); "
986  "L_{reco}/L_{truth} (leading); L_{reco}/L_{truth} (second)",
987  150,
988  0.,
989  1.5,
990  150,
991  0.,
992  1.5);
994 
996  tfs->make<TH2D>("h_MuonTrackStitching_MatchedCriteria_CriteriaTwoTracks_Angle",
997  "Good events: CriteriaTwoTracks vs. angle b/w leading and second muon track; "
998  "Criteria; Angle [#circ]",
999  7,
1000  0.75,
1001  4.25,
1002  100,
1003  0.,
1004  180.);
1006 
1007  //Stitching variables: bad events but leading plus second ok
1009  tfs->make<TH2D>(
1010  "h_MuonTrackStitching_FailedCriteriaButLeadingPlusSecondGood_TrackRes_Completeness",
1011  "Bad events but leading + second is good: L_{reco}/L_{truth} (leading) vs. Completeness "
1012  "(leading); L_{reco}/L_{truth} (leading); Completeness (leading)",
1013  150,
1014  0.,
1015  1.5,
1016  150,
1017  0.,
1018  1.5);
1020  "colz");
1021 
1023  tfs->make<TH2D>("h_MuonTrackStitching_FailedCriteriaButLeadingPlusSecondGood_Distance_Angle",
1024  "Bad events but leading + second is good: distance vs. angle b/w leading and "
1025  "second muon track; Distance [cm]; Angle [#circ]",
1026  100,
1027  0.,
1028  100.,
1029  100,
1030  0.,
1031  180.);
1033 
1035  tfs->make<TH2D>(
1036  "h_MuonTrackStitching_FailedCriteriaButLeadingPlusSecondGood_TrackResLeading_"
1037  "TrackResSecond",
1038  "Bad events but leading + second is good: L_{reco}/L_{truth} (leading) vs. "
1039  "L_{reco}/L_{truth} (second); L_{reco}/L_{truth} (leading); L_{reco}/L_{truth} (second)",
1040  150,
1041  0.,
1042  1.5,
1043  150,
1044  0.,
1045  1.5);
1047  ->SetOption("colz");
1048 
1049  //Stitching variables: bad events and leading plus second not ok
1051  tfs->make<TH2D>(
1052  "h_MuonTrackStitching_FailedCriteriaAndLeadingPlusSecondBad_TrackRes_Completeness",
1053  "Bad events and leading + second is bad: L_{reco}/L_{truth} (leading) vs. Completeness "
1054  "(leading); L_{reco}/L_{truth} (leading); Completeness (leading)",
1055  150,
1056  0.,
1057  1.5,
1058  150,
1059  0.,
1060  1.5);
1062  "colz");
1063 
1065  tfs->make<TH2D>("h_MuonTrackStitching_FailedCriteriaAndLeadingPlusSecondBad_Distance_Angle",
1066  "Bad events and leading + second is bad: distance vs. angle b/w leading and "
1067  "second muon track; Distance [cm]; Angle [#circ]",
1068  100,
1069  0.,
1070  100.,
1071  100,
1072  0.,
1073  180.);
1075 
1077  tfs->make<TH2D>(
1078  "h_MuonTrackStitching_FailedCriteriaAndLeadingPlusSecondBad_TrackResLeading_TrackResSecond",
1079  "Bad events and leading + second is bad: L_{reco}/L_{truth} (leading) vs. "
1080  "L_{reco}/L_{truth} (second); L_{reco}/L_{truth} (leading); L_{reco}/L_{truth} (second)",
1081  150,
1082  0.,
1083  1.5,
1084  150,
1085  0.,
1086  1.5);
1088  ->SetOption("colz");
1089  }
1090  //========================================================================
1091  void
1093  {
1094 
1095  doEfficiencies();
1096  }
1097  //========================================================================
1098  void
1100  {
1101  mf::LogInfo("MuonTrackingEff") << "begin run..." << std::endl;
1102  }
1103  //========================================================================
1104  void
1106  {
1107  if (event.isRealData()) return;
1108 
1109  bool isFiducial = false;
1110  processEff(event, isFiducial);
1111  }
1112  //========================================================================
1113  void
1115  {
1116 
1117  EventCounter++;
1118  simb::MCParticle* MCTruthMuonParticle = NULL;
1119 
1121  const sim::ParticleList& plist = pi_serv->ParticleList();
1122  simb::MCParticle* particle = 0;
1123 
1124  for (sim::ParticleList::const_iterator ipar = plist.begin(); ipar != plist.end(); ++ipar) {
1125  particle = ipar->second;
1126 
1127  if (particle->Mother() ==
1128  0) { //0=particle has no mother particle, 1=particle has a mother particle
1129  const TLorentzVector& positionStart = particle->Position(0);
1130  positionStart.GetXYZT(
1131  MCTruthMuonVertex); //MCTruthMuonVertex[0-2]=truth vertex, MCTruthMuonVertex[3]=t=0
1132  }
1133 
1134  if (particle->PdgCode() == fMuonPDGCode) { // Particle cannon muon
1135  MCTruthMuonParticle = particle;
1136  MCTruthMuonID = particle->TrackId();
1138  sqrt(pow(particle->Momentum().Px(), 2) + pow(particle->Momentum().Py(), 2) +
1139  pow(particle->Momentum().Pz(), 2));
1140 
1141  if (particle->Momentum().Pz() >= 0 && particle->Momentum().Px() >= 0) {
1143  (180.0 / 3.14159) * atan(particle->Momentum().Px() / particle->Momentum().Pz());
1144  }
1145  else if (particle->Momentum().Pz() < 0 && particle->Momentum().Px() >= 0) {
1147  180.0 + (180.0 / 3.14159) * atan(particle->Momentum().Px() / particle->Momentum().Pz());
1148  }
1149  else if (particle->Momentum().Pz() < 0 && particle->Momentum().Px() < 0) {
1151  180.0 + (180.0 / 3.14159) * atan(particle->Momentum().Px() / particle->Momentum().Pz());
1152  }
1153  else if (particle->Momentum().Pz() >= 0 && particle->Momentum().Px() < 0) {
1155  360.0 + (180.0 / 3.14159) * atan(particle->Momentum().Px() / particle->Momentum().Pz());
1156  }
1157 
1159  (180.0 / 3.14159) * asin(particle->Momentum().Py() / MCTruthMuonMomentum);
1160  }
1161  }
1162  double MCTruthLengthMuon = truthLength(MCTruthMuonParticle);
1163  h_TruthLength->Fill(MCTruthLengthMuon);
1164 
1165  //===================================================================
1166  //Saving denominator histograms
1167  //===================================================================
1168  isFiducial = insideFV(MCTruthMuonVertex);
1169  if (!isFiducial) return;
1170 
1171  //save events for Nucleon decay and particle cannon
1172  if (MCTruthMuonParticle) {
1175  h_SinThetaYZ_den->Fill(sin((3.14159 / 180.) * MCTruthMuonThetaYZ));
1178  sin((3.14159 / 180.) * MCTruthMuonThetaYZ));
1179  CountMCTruthMuon++;
1180  }
1181 
1182  //========================================================================
1183  // Reco stuff, once we have selected a MC Particle let's find out if there is a track associated
1184  //========================================================================
1185 
1186  int NMuonTracks = 0;
1187 
1188  art::Handle<std::vector<recob::Track>> TrackListHandle;
1189  if (!event.getByLabel(fTrackModuleLabel, TrackListHandle)) return;
1190  std::vector<art::Ptr<recob::Track>> TrackList;
1191  art::fill_ptr_vector(TrackList, TrackListHandle);
1192  int NRecoTracks = TrackList.size();
1193  art::FindManyP<recob::Hit> track_hits(TrackListHandle, event, fTrackModuleLabel);
1194  if (NRecoTracks == 0) {
1195  MF_LOG_DEBUG("MuonTrackingEff") << "There are no reco tracks... bye";
1196  std::cout << "There are no reco tracks! MCTruthMuonThetaXZ: " << std::endl;
1197 
1198  Criteria = 0.;
1201  sin((3.14159 / 180.) * MCTruthMuonThetaYZ));
1203  sin((3.14159 / 180.) * MCTruthMuonThetaYZ));
1204  h_Criteria_NRecoTrack_num->Fill(Criteria, static_cast<double>(NRecoTracks));
1205  h_Criteria_NMuonTrack_num->Fill(Criteria, static_cast<double>(NMuonTracks));
1207  return;
1208  }
1209  MF_LOG_DEBUG("MuonTrackingEff") << "Found this many reco tracks " << NRecoTracks;
1210 
1211  //std::cout << "NRecoTracks: " << NRecoTracks << std::endl;
1212 
1213  double PurityLeadingMuon = 0.;
1214  double CompletenessLeadingMuon = 0.;
1215  double RecoLengthLeadingMuon = 0.;
1216  art::Ptr<recob::Track> TrackLeadingMuon;
1217 
1218  double RecoLengthSecondMuon = 0.;
1219  double CompletenessSecondMuon = 0.;
1220  double PuritySecondMuon = 0.;
1221  art::Ptr<recob::Track> TrackSecondMuon;
1222 
1223  double TrackLengthMuonSum = 0.;
1224  double tmpTotalRecoEnergy = 0.;
1225 
1226  double MaxLengthNoRecoMuon = 0;
1227  int PDGCodeMaxLengthNoRecoMuon = 0;
1228 
1229  const simb::MCParticle* RecoMuonParticle = NULL;
1230 
1231  std::vector<art::Ptr<recob::Hit>> tmp_TrackHits = track_hits.at(0);
1232  std::vector<art::Ptr<recob::Hit>> AllHits;
1234  if (event.get(tmp_TrackHits[0].id(), HitHandle)) art::fill_ptr_vector(AllHits, HitHandle);
1235 
1236  auto const clockData =
1238 
1239  // Loop over reco tracks
1240  for (int i = 0; i < NRecoTracks; i++) {
1241  art::Ptr<recob::Track> track = TrackList[i];
1242  std::vector<art::Ptr<recob::Hit>> TrackHits = track_hits.at(i);
1243  double tmpPurity = 0.;
1244  double tmpCompleteness = 0.;
1245  const simb::MCParticle* particle;
1246 
1247  truthMatcher(
1248  clockData, AllHits, TrackHits, particle, tmpPurity, tmpCompleteness, tmpTotalRecoEnergy);
1249 
1250  if (!particle) {
1251  std::cout << "ERROR: Truth matcher didn't find a particle!" << std::endl;
1252  continue;
1253  }
1254 
1255  if (track->Length() > MaxLengthNoRecoMuon && particle->PdgCode() != fMuonPDGCode &&
1256  particle->TrackId() != MCTruthMuonID) {
1257  MaxLengthNoRecoMuon = track->Length();
1258  PDGCodeMaxLengthNoRecoMuon = particle->PdgCode();
1259  }
1260  //some muon tracks have Completeness=0 and Purity=0 but are still considered as muon tracks in function truthmatcher. Getting rid of these tracks by asking tmpCompleteness > 0 && tmpPurity > 0
1261  if ((particle->PdgCode() == fMuonPDGCode) && (particle->TrackId() == MCTruthMuonID) &&
1262  tmpCompleteness > 0 && tmpPurity > 0) {
1263 
1264  NMuonTracks++;
1265  TrackLengthMuonSum += track->Length();
1266 
1267  if (NMuonTracks == 1) {
1268  CompletenessLeadingMuon = tmpCompleteness;
1269  PurityLeadingMuon = tmpPurity;
1270  RecoLengthLeadingMuon = track->Length();
1271  TrackLeadingMuon = track;
1272 
1273  RecoMuonParticle = particle;
1274  }
1275 
1276  if (NMuonTracks >= 2 && tmpCompleteness > CompletenessLeadingMuon) {
1277 
1278  CompletenessSecondMuon = CompletenessLeadingMuon;
1279  PuritySecondMuon = PurityLeadingMuon;
1280  RecoLengthSecondMuon = RecoLengthLeadingMuon;
1281  TrackSecondMuon = TrackLeadingMuon;
1282 
1283  CompletenessLeadingMuon = tmpCompleteness;
1284  PurityLeadingMuon = tmpPurity;
1285  RecoLengthLeadingMuon = track->Length();
1286  TrackLeadingMuon = track;
1287 
1288  RecoMuonParticle = particle;
1289  }
1290 
1291  else if (NMuonTracks >= 2 && tmpCompleteness < CompletenessLeadingMuon &&
1292  tmpCompleteness > CompletenessSecondMuon) {
1293  CompletenessSecondMuon = tmpCompleteness;
1294  PuritySecondMuon = tmpPurity;
1295  RecoLengthSecondMuon = track->Length();
1296  TrackSecondMuon = track;
1297  }
1298  }
1299  }
1300 
1301  h_TotalRecoEnergy->Fill(tmpTotalRecoEnergy);
1302 
1303  //Muon events
1304  //=========================================================
1305 
1306  if (RecoMuonParticle && MCTruthMuonParticle) {
1307  CountRecoMuon++;
1308  h_Purity->Fill(PurityLeadingMuon);
1309  h_Completeness->Fill(CompletenessLeadingMuon);
1310  h_TrackRes->Fill(RecoLengthLeadingMuon / MCTruthLengthMuon);
1311 
1312  std::cout << "TrackLeadingMuon->Vertex().X(): " << TrackLeadingMuon->Vertex().X()
1313  << std::endl;
1314  std::cout << "MCTruthMuonParticle->Vz(): " << MCTruthMuonParticle->Vz() << std::endl;
1315  std::cout << " " << std::endl;
1316 
1317  double DistanceBetweenTruthAndRecoTrack;
1318  double AngleBetweenTruthAndRecoTrack;
1320  TrackLeadingMuon,
1321  DistanceBetweenTruthAndRecoTrack,
1322  AngleBetweenTruthAndRecoTrack);
1323 
1324  h_VertexRes->Fill(DistanceBetweenTruthAndRecoTrack);
1325  h_DirectionRes->Fill(AngleBetweenTruthAndRecoTrack);
1326 
1327  h_MuonTrackStitching_TrackRes_Completeness->Fill(RecoLengthLeadingMuon / MCTruthLengthMuon,
1328  CompletenessLeadingMuon);
1329  if (NMuonTracks >= 2) {
1330  double DistanceBetweenTracks;
1331  double AngleBetweenTracks;
1332  double CriteriaTwoTracks;
1333 
1334  FuncDistanceAndAngleBetweenTracks(TrackSecondMuon,
1335  TrackLeadingMuon,
1336  DistanceBetweenTracks,
1337  AngleBetweenTracks,
1338  CriteriaTwoTracks);
1339 
1341  RecoLengthLeadingMuon / MCTruthLengthMuon, RecoLengthSecondMuon / MCTruthLengthMuon);
1342  h_MuonTrackStitching_Distance_Angle->Fill(DistanceBetweenTracks, AngleBetweenTracks);
1344  RecoLengthSecondMuon / MCTruthLengthMuon, AngleBetweenTracks);
1345  h_MuonTrackStitching_CompletenessSecondMuon_Angle->Fill(CompletenessSecondMuon,
1346  AngleBetweenTracks);
1347  h_MuonTrackStitching_CriteriaTwoTracks_Angle->Fill(CriteriaTwoTracks, AngleBetweenTracks);
1348  }
1349 
1350  //Completeness
1351  if (CompletenessLeadingMuon < 0.5) {
1353  Criteria = 2.;
1355  sin((3.14159 / 180.) * MCTruthMuonThetaYZ));
1356  h_Criteria_NRecoTrack_num->Fill(Criteria, static_cast<double>(NRecoTracks));
1357  h_Criteria_NMuonTrack_num->Fill(Criteria, static_cast<double>(NMuonTracks));
1358  }
1359  //Purity
1360  if (PurityLeadingMuon < 0.5) {
1361  CountPurity++;
1362  Criteria = 3.;
1364  sin((3.14159 / 180.) * MCTruthMuonThetaYZ));
1365  h_Criteria_NRecoTrack_num->Fill(Criteria, static_cast<double>(NRecoTracks));
1366  h_Criteria_NMuonTrack_num->Fill(Criteria, static_cast<double>(NMuonTracks));
1367  }
1368  //Track too short
1369  if (RecoLengthLeadingMuon / MCTruthLengthMuon < 0.75) {
1371  Criteria = 4.;
1373  sin((3.14159 / 180.) * MCTruthMuonThetaYZ));
1374  h_Criteria_NRecoTrack_num->Fill(Criteria, static_cast<double>(NRecoTracks));
1375  h_Criteria_NMuonTrack_num->Fill(Criteria, static_cast<double>(NMuonTracks));
1376  }
1377  //Track too long
1378  if (RecoLengthLeadingMuon / MCTruthLengthMuon > 1.25) {
1380  Criteria = 5.;
1382  sin((3.14159 / 180.) * MCTruthMuonThetaYZ));
1383  h_Criteria_NRecoTrack_num->Fill(Criteria, static_cast<double>(NRecoTracks));
1384  h_Criteria_NMuonTrack_num->Fill(Criteria, static_cast<double>(NMuonTracks));
1385 
1386  /*std::cout << "Track too long!" << std::endl;
1387  art::ServiceHandle<cheat::BackTracker const> bt2;
1388  const sim::ParticleList& plist2 = bt2->ParticleList();
1389  simb::MCParticle *particle2=0;
1390 
1391  std::cout << "EventCounter: " << EventCounter << std::endl;
1392  for( sim::ParticleList::const_iterator ipar2 = plist2.begin(); ipar2!=plist2.end(); ++ipar2)
1393  {
1394  particle2 = ipar2->second;
1395  std::cout << "particle2->TrackId(): " << particle2->TrackId() << std::endl;
1396  std::cout << "particle2->PdgCode(): " << particle2->PdgCode() << std::endl;
1397  std::cout << "particle2->Momentum().P(): " << particle2->Momentum().P() << std::endl;
1398  std::cout << "tuthLength(particle2): " << truthLength(particle2) << std::endl;
1399  std::cout << "particle2->Position(): (x,y,z): " << particle2->Vx() << "\t" << particle2->Vy() << "\t" << particle2->Vz() << std::endl;
1400  std::cout << "particle2->Momentum(): (Px,Py,Pz): " << particle2->Momentum().Px() << "\t" << particle2->Momentum().Py() << "\t" << particle2->Momentum().Pz() << std::endl;
1401  std::cout << "particle2->Position().T(): " << particle2->Position().T() << std::endl;
1402  std::cout << "" << std::endl;
1403  }*/
1404  }
1405  //Reco failed at least one of the above criteria
1406  if (CompletenessLeadingMuon < 0.5 || PurityLeadingMuon < 0.5 ||
1407  RecoLengthLeadingMuon / MCTruthLengthMuon < 0.75 ||
1408  RecoLengthLeadingMuon / MCTruthLengthMuon > 1.25) {
1412  sin((3.14159 / 180.) * MCTruthMuonThetaYZ));
1414  RecoLengthLeadingMuon / MCTruthLengthMuon, CompletenessLeadingMuon);
1415 
1416  if (NMuonTracks == 1) { BadEvents1MuonTrack++; }
1417  if (NMuonTracks == 2) { BadEvents2MuonTrack++; }
1418  if (NMuonTracks == 3) { BadEvents3MuonTrack++; }
1419  if (NMuonTracks == 4) { BadEvents4OrMoreMuonTrack++; }
1420 
1421  if (NMuonTracks >= 2) {
1422  double AngleBetweenTracks;
1423  double DistanceBetweenTracks;
1424  double CriteriaTwoTracks;
1425  FuncDistanceAndAngleBetweenTracks(TrackSecondMuon,
1426  TrackLeadingMuon,
1427  DistanceBetweenTracks,
1428  AngleBetweenTracks,
1429  CriteriaTwoTracks);
1430 
1431  if (AngleBetweenTracks > 160.) {
1432  std::cout << "EventCounter: " << EventCounter << std::endl;
1433  std::cout << "Angle b/w tracks: " << AngleBetweenTracks << std::endl;
1434  std::cout << "CriteriaTwoTracks: " << CriteriaTwoTracks << std::endl;
1435  std::cout << "CompletenessLeadingMuon: " << CompletenessLeadingMuon << std::endl;
1436  std::cout << "CompletenessSecondMuon: " << CompletenessSecondMuon << std::endl;
1437  std::cout << "PurityLeadingMuon: " << PurityLeadingMuon << std::endl;
1438  std::cout << "PuritySecondMuon: " << PuritySecondMuon << std::endl;
1439  std::cout << "TrackLeadingMuon: " << RecoLengthLeadingMuon / MCTruthLengthMuon
1440  << std::endl;
1441  std::cout << "TrackResSecondMuon: " << RecoLengthSecondMuon / MCTruthLengthMuon
1442  << std::endl;
1443  std::cout << "TrackID leading: " << TrackLeadingMuon->ID() << std::endl;
1444  std::cout << "TrackID second: " << TrackSecondMuon->ID() << std::endl;
1445  }
1446 
1447  h_MuonTrackStitching_FailedCriteria_Distance_Angle->Fill(DistanceBetweenTracks,
1448  AngleBetweenTracks);
1450  RecoLengthLeadingMuon / MCTruthLengthMuon, RecoLengthSecondMuon / MCTruthLengthMuon);
1452  CompletenessLeadingMuon, CompletenessSecondMuon);
1454  RecoLengthSecondMuon / MCTruthLengthMuon, AngleBetweenTracks);
1456  CompletenessSecondMuon, AngleBetweenTracks);
1458  AngleBetweenTracks);
1459  if ((CompletenessLeadingMuon + CompletenessSecondMuon) >= 0.5 &&
1460  PurityLeadingMuon >= 0.5 && PuritySecondMuon >= 0.5 &&
1461  (RecoLengthLeadingMuon + RecoLengthSecondMuon) / MCTruthLengthMuon >= 0.75 &&
1462  (RecoLengthLeadingMuon + RecoLengthSecondMuon) / MCTruthLengthMuon <= 1.25) {
1466  MCTruthMuonThetaXZ, sin((3.14159 / 180.) * MCTruthMuonThetaYZ));
1467 
1469  RecoLengthLeadingMuon / MCTruthLengthMuon, CompletenessLeadingMuon);
1471  ->Fill(RecoLengthLeadingMuon / MCTruthLengthMuon,
1472  RecoLengthSecondMuon / MCTruthLengthMuon);
1474  DistanceBetweenTracks, AngleBetweenTracks);
1475  }
1476  if ((CompletenessLeadingMuon + CompletenessSecondMuon) < 0.5 || PurityLeadingMuon < 0.5 ||
1477  PuritySecondMuon < 0.5 ||
1478  (RecoLengthLeadingMuon + RecoLengthSecondMuon) / MCTruthLengthMuon < 0.75 ||
1479  (RecoLengthLeadingMuon + RecoLengthSecondMuon) / MCTruthLengthMuon > 1.25) {
1481 
1483  RecoLengthLeadingMuon / MCTruthLengthMuon, CompletenessLeadingMuon);
1485  ->Fill(RecoLengthLeadingMuon / MCTruthLengthMuon,
1486  RecoLengthSecondMuon / MCTruthLengthMuon);
1488  DistanceBetweenTracks, AngleBetweenTracks);
1489  }
1490  if ((CompletenessLeadingMuon + CompletenessSecondMuon) < 0.5) {
1492  }
1493  if (PurityLeadingMuon < 0.5 || PuritySecondMuon < 0.5) {
1495  }
1496  if ((RecoLengthLeadingMuon + RecoLengthSecondMuon) / MCTruthLengthMuon < 0.75) {
1498  }
1499  if ((RecoLengthLeadingMuon + RecoLengthSecondMuon) / MCTruthLengthMuon > 1.25) {
1501  }
1502  }
1503  else if (NMuonTracks == 1) {
1505  if (CompletenessLeadingMuon < 0.5) {
1507  }
1508  if (PurityLeadingMuon < 0.5) { CountBadLeadingMuonTrackAndOnlyOneMuonTrackPurity++; }
1509  if (RecoLengthLeadingMuon / MCTruthLengthMuon < 0.75) {
1511  }
1512  if (RecoLengthLeadingMuon / MCTruthLengthMuon > 1.25) {
1514  }
1515  }
1516  }
1517  //Reco succesful according to the above criteria
1518  if (CompletenessLeadingMuon >= 0.5 && PurityLeadingMuon >= 0.5 &&
1519  RecoLengthLeadingMuon / MCTruthLengthMuon >= 0.75 &&
1520  RecoLengthLeadingMuon / MCTruthLengthMuon <= 1.25) {
1522  Criteria = 6.;
1525  h_SinThetaYZ_num->Fill(sin((3.14159 / 180.) * MCTruthMuonThetaYZ));
1528  sin((3.14159 / 180.) * MCTruthMuonThetaYZ));
1529  h_Criteria_NRecoTrack_num->Fill(Criteria, static_cast<double>(NRecoTracks));
1530  h_Criteria_NMuonTrack_num->Fill(Criteria, static_cast<double>(NMuonTracks));
1531 
1533  RecoLengthLeadingMuon / MCTruthLengthMuon, CompletenessLeadingMuon);
1534 
1535  if (NMuonTracks == 1) { GoodEvents1MuonTrack++; }
1536  if (NMuonTracks == 2) { GoodEvents2MuonTrack++; }
1537  if (NMuonTracks == 3) { GoodEvents3MuonTrack++; }
1538  if (NMuonTracks == 4) { GoodEvents4OrMoreMuonTrack++; }
1539 
1540  if (NMuonTracks >= 2) {
1542  RecoLengthLeadingMuon / MCTruthLengthMuon, RecoLengthSecondMuon / MCTruthLengthMuon);
1543 
1544  double AngleBetweenTracks;
1545  double DistanceBetweenTracks;
1546  double CriteriaTwoTracks;
1547  FuncDistanceAndAngleBetweenTracks(TrackSecondMuon,
1548  TrackLeadingMuon,
1549  DistanceBetweenTracks,
1550  AngleBetweenTracks,
1551  CriteriaTwoTracks);
1552  h_MuonTrackStitching_MatchedCriteria_Distance_Angle->Fill(DistanceBetweenTracks,
1553  AngleBetweenTracks);
1555  AngleBetweenTracks);
1556  }
1557  }
1558  }
1559  //No muon track
1560  if (!RecoMuonParticle && MCTruthMuonParticle) {
1563  Criteria = 1.;
1565  sin((3.14159 / 180.) * MCTruthMuonThetaYZ));
1568  sin((3.14159 / 180.) * MCTruthMuonThetaYZ));
1569  h_Criteria_NRecoTrack_num->Fill(Criteria, static_cast<double>(NRecoTracks));
1570  h_Criteria_NMuonTrack_num->Fill(Criteria, static_cast<double>(NMuonTracks));
1571  h_NoMuonTrack_MaxTrackLength_PDGCode->Fill(MaxLengthNoRecoMuon,
1572  static_cast<double>(PDGCodeMaxLengthNoRecoMuon));
1573  }
1574  }
1575  //========================================================================
1576  void
1579  std::vector<art::Ptr<recob::Hit>> track_hits,
1580  const simb::MCParticle*& MCparticle,
1581  double& Purity,
1582  double& Completeness,
1583  double& TotalRecoEnergy)
1584  {
1587  std::map<int, double> trkID_E; // map that connects TrackID and energy for
1588  // each hit <trackID, energy>
1589  for (size_t j = 0; j < track_hits.size(); ++j) {
1590  art::Ptr<recob::Hit> hit = track_hits[j];
1591  std::vector<sim::TrackIDE> TrackIDs =
1592  bt_serv->HitToTrackIDEs(clockData,
1593  hit); // TrackIDE contains TrackID, energy and energyFrac. A hit can
1594  // have several TrackIDs (so this hit is associated with multiple
1595  // MC truth track IDs (EM shower IDs are negative). If a hit ahs
1596  // multiple trackIDs, "energyFrac" contains the fraction of the
1597  // energy of for each ID compared to the total energy of the hit.
1598  // "energy" contains only the energy associated with the specific
1599  // ID in that case. This requires MC truth info!
1600  for (size_t k = 0; k < TrackIDs.size(); k++) {
1601  trkID_E[TrackIDs[k].trackID] +=
1602  TrackIDs[k].energy; // sum up the energy for each TrackID and store
1603  // <TrackID, energy> in "TrkID_E"
1604  }
1605  }
1606 
1607  double E_em = 0.0;
1608  double max_E = -999.0;
1609  double TotalEnergyTrack = 0.0;
1610  int TrackID = -999;
1611  double PartialEnergyTrackID =
1612  0.0; // amount of energy deposited by the particle that deposited more
1613  // energy... tomato potato... blabla
1614  //! if the collection of hits have more than one particle associate save the
1615  //! particle w/ the highest energy deposition since we are looking for
1616  //! muons/pions/protons this should be enough
1617  if (!trkID_E.size()) {
1618  MCparticle = 0;
1619  return; // Ghost track???
1620  }
1621  for (std::map<int, double>::iterator ii = trkID_E.begin(); ii != trkID_E.end();
1622  ++ii) { // trkID_E contains the trackID (first) and corresponding
1623  // energy (second) for a specific track, summed up over all
1624  // events. here looping over all trekID_E's
1625  TotalEnergyTrack += ii->second; // and summing up the energy of all hits
1626  // in the track (TotalEnergyTrack)
1627  if ((ii->second) > max_E) { // looking for the trakID with the highest energy in the
1628  // track. this is PartialEnergyTrackID and max_E then
1629  PartialEnergyTrackID = ii->second;
1630  max_E = ii->second;
1631  TrackID = ii->first; // saving trackID of the ID with the highest energy
1632  // contribution in the track to assign it to
1633  // MCparticle later
1634  if (TrackID < 0) E_em += ii->second; // IDs of em shower particles are negative
1635  }
1636  }
1637  MCparticle = pi_serv->TrackIdToParticle_P(TrackID);
1638  // In the current simulation, we do not save EM Shower daughters in GEANT.
1639  // But we do save the energy deposition in TrackIDEs. If the energy
1640  // deposition is from a particle that is the daughter of an EM particle, the
1641  // negative of the parent track ID is saved in TrackIDE for the daughter
1642  // particle we don't want to track gammas or any other EM activity
1643  if (TrackID < 0) { return; }
1644 
1645  // Purity = (PartialEnergyTrackID+E_em)/TotalEnergyTrack;
1646  Purity = PartialEnergyTrackID / TotalEnergyTrack;
1647 
1648  // completeness
1649  TotalRecoEnergy = 0;
1650  for (size_t k = 0; k < AllHits.size();
1651  ++k) { // loop over all hits (all hits in all tracks of the event, not
1652  // only the hits in the track we were looking at before)
1653  art::Ptr<recob::Hit> hit = AllHits[k];
1654  std::vector<sim::TrackIDE> TrackIDs = bt_serv->HitToTrackIDEs(clockData, hit);
1655  for (size_t l = 0; l < TrackIDs.size(); ++l) { // and over all track IDs of the hits
1656  if (TrackIDs[l].trackID == TrackID)
1657  TotalRecoEnergy += TrackIDs[l].energy; // and sum up the energy fraction of all hits
1658  // that correspond ot the saved trackID
1659  }
1660  }
1661  Completeness = PartialEnergyTrackID / TotalRecoEnergy;
1662  }
1663 
1664  void
1666  art::Ptr<recob::Track> Track2,
1667  double& TempDistanceBetweenTracks,
1668  double& TempAngleBetweenTracks,
1669  double& TempCriteriaTwoTracks)
1670  {
1671 
1672  TempDistanceBetweenTracks = sqrt(pow(Track1->End().X() - Track2->Vertex().X(), 2) +
1673  pow(Track1->End().Y() - Track2->Vertex().Y(), 2) +
1674  pow(Track1->End().Z() - Track2->Vertex().Z(), 2));
1675  TempAngleBetweenTracks = (180.0 / 3.14159) * Track1->EndDirection<TVector3>().Angle(
1676  Track2->VertexDirection<TVector3>());
1677  TempCriteriaTwoTracks = 1.;
1678 
1679  if (TempDistanceBetweenTracks > sqrt(pow(Track1->End().X() - Track2->End().X(), 2) +
1680  pow(Track1->End().Y() - Track2->End().Y(), 2) +
1681  pow(Track1->End().Z() - Track2->End().Z(), 2))) {
1682  TempDistanceBetweenTracks = sqrt(pow(Track1->End().X() - Track2->End().X(), 2) +
1683  pow(Track1->End().Y() - Track2->End().Y(), 2) +
1684  pow(Track1->End().Z() - Track2->End().Z(), 2));
1685  TempAngleBetweenTracks = 180. - (180.0 / 3.14159) * Track1->EndDirection<TVector3>().Angle(
1686  Track2->EndDirection<TVector3>());
1687  TempCriteriaTwoTracks = 2.;
1688  }
1689 
1690  if (TempDistanceBetweenTracks > sqrt(pow(Track1->Vertex().X() - Track2->End().X(), 2) +
1691  pow(Track1->Vertex().Y() - Track2->End().Y(), 2) +
1692  pow(Track1->Vertex().Z() - Track2->End().Z(), 2))) {
1693  TempDistanceBetweenTracks = sqrt(pow(Track1->Vertex().X() - Track2->End().X(), 2) +
1694  pow(Track1->Vertex().Y() - Track2->End().Y(), 2) +
1695  pow(Track1->Vertex().Z() - Track2->End().Z(), 2));
1696  TempAngleBetweenTracks = (180.0 / 3.14159) * Track1->VertexDirection<TVector3>().Angle(
1697  Track2->EndDirection<TVector3>());
1698  TempCriteriaTwoTracks = 3.;
1699  }
1700 
1701  if (TempDistanceBetweenTracks > sqrt(pow(Track1->Vertex().X() - Track2->Vertex().X(), 2) +
1702  pow(Track1->Vertex().Y() - Track2->Vertex().Y(), 2) +
1703  pow(Track1->Vertex().Z() - Track2->Vertex().Z(), 2))) {
1704  TempDistanceBetweenTracks = sqrt(pow(Track1->Vertex().X() - Track2->Vertex().X(), 2) +
1705  pow(Track1->Vertex().Y() - Track2->Vertex().Y(), 2) +
1706  pow(Track1->Vertex().Z() - Track2->Vertex().Z(), 2));
1707  TempAngleBetweenTracks = 180. - (180.0 / 3.14159) * Track1->VertexDirection<TVector3>().Angle(
1708  Track2->VertexDirection<TVector3>());
1709  TempCriteriaTwoTracks = 4.;
1710  }
1711  }
1712 
1713  void
1715  const simb::MCParticle*& MCparticle,
1717  double& TempDistanceBetweenTruthAndRecoTrack,
1718  double& TempAngleBeetweenTruthAndRecoTrack)
1719  {
1720  TempDistanceBetweenTruthAndRecoTrack = sqrt(pow(Track->Vertex().X() - MCparticle->Vx(), 2) +
1721  pow(Track->Vertex().Y() - MCparticle->Vy(), 2) +
1722  pow(Track->Vertex().Z() - MCparticle->Vz(), 2));
1723 
1724  TempAngleBeetweenTruthAndRecoTrack = 0;
1725  }
1726 
1727  //========================================================================
1728  double
1730  {
1731  //calculate the truth length considering only the part that is inside the TPC
1732  //Base on a piece of code from dune/TrackingAna/TrackingEfficiency_module.cc
1733 
1734  if (!MCparticle) return -999.0;
1735  int numberTrajectoryPoints = MCparticle->NumberTrajectoryPoints();
1736  std::vector<double> TPCLengthHits(numberTrajectoryPoints, 0);
1737  int FirstHit = 0, LastHit = 0;
1738  double TPCLength = 0.0;
1739  bool BeenInVolume = false;
1740 
1741  for (unsigned int MCHit = 0; MCHit < TPCLengthHits.size(); ++MCHit) {
1742  const TLorentzVector& tmpPosition = MCparticle->Position(MCHit);
1743  double const tmpPosArray[] = {tmpPosition[0], tmpPosition[1], tmpPosition[2]};
1744  if (MCHit != 0)
1745  TPCLengthHits[MCHit] = sqrt(pow((MCparticle->Vx(MCHit - 1) - MCparticle->Vx(MCHit)), 2) +
1746  pow((MCparticle->Vy(MCHit - 1) - MCparticle->Vy(MCHit)), 2) +
1747  pow((MCparticle->Vz(MCHit - 1) - MCparticle->Vz(MCHit)), 2));
1748  geo::TPCID tpcid = geom->FindTPCAtPosition(tmpPosArray);
1749  if (tpcid.isValid) {
1750  LastHit = MCHit;
1751  if (!BeenInVolume) {
1752  BeenInVolume = true;
1753  FirstHit = MCHit;
1754  }
1755  }
1756  }
1757  for (int Hit = FirstHit + 1; Hit <= LastHit; ++Hit)
1758  TPCLength += TPCLengthHits[Hit];
1759  return TPCLength;
1760  }
1761  //========================================================================
1762  bool
1764  {
1765 
1766  double x = vertex[0];
1767  double y = vertex[1];
1768  double z = vertex[2];
1769 
1770  if (x > fFidVolXmin && x < fFidVolXmax && y > fFidVolYmin && y < fFidVolYmax &&
1771  z > fFidVolZmin && z < fFidVolZmax)
1772  return true;
1773  else
1774  return false;
1775  }
1776  //========================================================================
1777  void
1779  {
1780  std::cout << std::endl;
1781 
1782  std::cout << "EventCounter: "
1783  << "\t" << EventCounter << std::endl;
1784 
1785  std::cout << "CountMCTruthMuon: "
1786  << "\t" << CountMCTruthMuon << " = "
1787  << 100 * static_cast<double>(CountMCTruthMuon) / static_cast<double>(EventCounter)
1788  << "%" << std::endl;
1789 
1790  std::cout << "CountGoodLeadingMuonTrack (=good events): "
1791  << "\t" << CountGoodLeadingMuonTrack << "/" << CountMCTruthMuon << " = "
1792  << 100 * static_cast<double>(CountGoodLeadingMuonTrack) /
1793  static_cast<double>(CountMCTruthMuon)
1794  << "%" << std::endl;
1795 
1796  std::cout << "CountBadLeadingMuonTrack+CountNoRecoTracks+CountNoMuonTracks (=bad events): "
1798  << 100 *
1799  static_cast<double>(CountBadLeadingMuonTrack + CountNoRecoTracks +
1801  static_cast<double>(CountMCTruthMuon)
1802  << "%" << std::endl;
1803 
1804  std::cout << "CountNoRecoTracks+CountNoMuonTracks: "
1805  << "\t" << CountNoRecoTracks + CountNoMuonTracks << " = "
1806  << 100 * static_cast<double>(CountNoRecoTracks + CountNoMuonTracks) /
1807  static_cast<double>(CountMCTruthMuon)
1808  << "%" << std::endl;
1809 
1810  std::cout << "CountTrackLengthTooShort: "
1811  << "\t" << CountTrackLengthTooShort << " = "
1812  << 100 * static_cast<double>(CountTrackLengthTooShort) /
1813  static_cast<double>(CountMCTruthMuon)
1814  << "%" << std::endl;
1815 
1816  std::cout << "CountCompleteness: "
1817  << "\t" << CountCompleteness << " = "
1818  << 100 * static_cast<double>(CountCompleteness) /
1819  static_cast<double>(CountMCTruthMuon)
1820  << "%" << std::endl;
1821 
1822  std::cout << "CountTrackLengthTooLong: "
1823  << "\t" << CountTrackLengthTooLong << " = "
1824  << 100 * static_cast<double>(CountTrackLengthTooLong) /
1825  static_cast<double>(CountMCTruthMuon)
1826  << "%" << std::endl;
1827 
1828  std::cout << "CountPurity: "
1829  << "\t" << CountPurity << " = "
1830  << 100 * static_cast<double>(CountPurity) / static_cast<double>(CountMCTruthMuon)
1831  << "%" << std::endl;
1832 
1833  std::cout << std::endl;
1834 
1835  std::cout << "GoodLeadingMuonTrack+CountBadLeadingMuonTrackButLeadingPlusSecondGood (=good "
1836  "events after stitching): "
1837  << "\t"
1839  << CountMCTruthMuon << " = "
1840  << 100 *
1841  static_cast<double>(CountGoodLeadingMuonTrack +
1843  static_cast<double>(CountMCTruthMuon)
1844  << "%" << std::endl;
1845 
1846  std::cout << "CountBadLeadingMuonTrack+CountNoRecoTracks+CountNoMuonTracks-"
1847  "CountBadLeadingMuonTrackButLeadingPlusSecondGood (=bad events after stitching) : "
1848  << "\t"
1851  << " = "
1852  << 100 *
1853  static_cast<double>(CountBadLeadingMuonTrack + CountNoRecoTracks +
1856  static_cast<double>(CountMCTruthMuon)
1857  << "%" << std::endl;
1858 
1859  std::cout << std::endl;
1860 
1861  std::cout << "CountBadLeadingMuonTrackButLeadingPlusSecondGood: "
1863  << 100 * static_cast<double>(CountBadLeadingMuonTrackButLeadingPlusSecondGood) /
1864  static_cast<double>(CountMCTruthMuon)
1865  << "%" << std::endl;
1866 
1867  std::cout << "CountBadLeadingMuonTrackAndOnlyOneMuonTrack: "
1869  << 100 * static_cast<double>(CountBadLeadingMuonTrackAndOnlyOneMuonTrack) /
1870  static_cast<double>(CountMCTruthMuon)
1871  << "%" << std::endl;
1872 
1873  std::cout << "CountBadLeadingMuonTrackAndLeadingPlusSecondBad: "
1875  << 100 * static_cast<double>(CountBadLeadingMuonTrackAndLeadingPlusSecondBad) /
1876  static_cast<double>(CountMCTruthMuon)
1877  << "%" << std::endl;
1878 
1879  std::cout << "CountNoRecoTracks: "
1880  << "\t" << CountNoRecoTracks << " = "
1881  << 100 * static_cast<double>(CountNoRecoTracks) /
1882  static_cast<double>(CountMCTruthMuon)
1883  << "%" << std::endl;
1884 
1885  std::cout << "CountNoMuonTracks: "
1886  << "\t" << CountNoMuonTracks << " = "
1887  << 100 * static_cast<double>(CountNoMuonTracks) /
1888  static_cast<double>(CountMCTruthMuon)
1889  << "%" << std::endl;
1890 
1891  std::cout << std::endl;
1892 
1893  std::cout << "CountBadLeadingMuonTrackAndOnlyOneMuonTrackCompleteness: "
1895  std::cout << "CountBadLeadingMuonTrackAndOnlyOneMuonTrackPurity: "
1897  std::cout << "CountBadLeadingMuonTrackAndOnlyOneMuonTrackTrackTooShort: "
1899  std::cout << "CountBadLeadingMuonTrackAndOnlyOneMuonTrackTrackTooLong: "
1901 
1902  std::cout << std::endl;
1903 
1904  std::cout << "CountBadLeadingMuonTrackAndLeadingPlusSecondBadCompleteness: "
1906  << 100 *
1907  static_cast<double>(
1909  static_cast<double>(CountMCTruthMuon)
1910  << "%" << std::endl;
1911 
1912  std::cout << "CountBadLeadingMuonTrackAndLeadingPlusSecondBadPurity: "
1914  << 100 * static_cast<double>(CountBadLeadingMuonTrackAndLeadingPlusSecondBadPurity) /
1915  static_cast<double>(CountMCTruthMuon)
1916  << "%" << std::endl;
1917 
1918  std::cout << "CountBadLeadingMuonTrackAndLeadingPlusSecondBadTrackTooShort: "
1920  << 100 *
1921  static_cast<double>(
1923  static_cast<double>(CountMCTruthMuon)
1924  << "%" << std::endl;
1925 
1926  std::cout << "CountBadLeadingMuonTrackAndLeadingPlusSecondBadTrackTooLong: "
1928  << 100 *
1929  static_cast<double>(
1931  static_cast<double>(CountMCTruthMuon)
1932  << "%" << std::endl;
1933 
1934  std::cout << std::endl;
1935 
1936  std::cout << "GoodEvents1MuonTrack: " << GoodEvents1MuonTrack << std::endl;
1937  std::cout << "GoodEvents2MuonTrack: " << GoodEvents2MuonTrack << std::endl;
1938  std::cout << "GoodEvents3MuonTrack: " << GoodEvents3MuonTrack << std::endl;
1939  std::cout << "GoodEvents4OrMoreMuonTrack: " << GoodEvents4OrMoreMuonTrack << std::endl;
1940 
1941  std::cout << "BadEvents0MuonTrack: " << BadEvents0MuonTrack << std::endl;
1942  std::cout << "BadEvents1MuonTrack: " << BadEvents1MuonTrack << std::endl;
1943  std::cout << "BadEvents2MuonTrack: " << BadEvents2MuonTrack << std::endl;
1944  std::cout << "BadEvents3MuonTrack: " << BadEvents3MuonTrack << std::endl;
1945  std::cout << "BadEvents4OrMoreMuonTrack: " << BadEvents4OrMoreMuonTrack << std::endl;
1946 
1948 
1949  for (int i = 0; i < (NCriteriaBins + 1) / 2; ++i) {
1950  for (int j = 0; j < (NRecoTracksBins + 1) / 2; ++j) {
1951  if (i == 0) {
1952  h_Criteria_NRecoTrack_den->SetBinContent(1 + 2 * i, 1 + 2 * j, CountNoRecoTracks);
1953  h_Criteria_NMuonTrack_den->SetBinContent(1 + 2 * i, 1 + 2 * j, CountNoRecoTracks);
1954  }
1955  else if (i == 1) {
1956  h_Criteria_NRecoTrack_den->SetBinContent(1 + 2 * i, 1 + 2 * j, CountNoMuonTracks);
1957  h_Criteria_NMuonTrack_den->SetBinContent(1 + 2 * i, 1 + 2 * j, CountNoMuonTracks);
1958  }
1959  else if (i == 2) {
1960  h_Criteria_NRecoTrack_den->SetBinContent(1 + 2 * i, 1 + 2 * j, CountCompleteness);
1961  h_Criteria_NMuonTrack_den->SetBinContent(1 + 2 * i, 1 + 2 * j, CountCompleteness);
1962  }
1963  else if (i == 3) {
1964  h_Criteria_NRecoTrack_den->SetBinContent(1 + 2 * i, 1 + 2 * j, CountPurity);
1965  h_Criteria_NMuonTrack_den->SetBinContent(1 + 2 * i, 1 + 2 * j, CountPurity);
1966  }
1967  else if (i == 4) {
1968  h_Criteria_NRecoTrack_den->SetBinContent(1 + 2 * i, 1 + 2 * j, CountTrackLengthTooShort);
1969  h_Criteria_NMuonTrack_den->SetBinContent(1 + 2 * i, 1 + 2 * j, CountTrackLengthTooShort);
1970  }
1971  else if (i == 5) {
1972  h_Criteria_NRecoTrack_den->SetBinContent(1 + 2 * i, 1 + 2 * j, CountTrackLengthTooLong);
1973  h_Criteria_NMuonTrack_den->SetBinContent(1 + 2 * i, 1 + 2 * j, CountTrackLengthTooLong);
1974  }
1975  else if (i == 6) {
1976  h_Criteria_NRecoTrack_den->SetBinContent(1 + 2 * i, 1 + 2 * j, CountRecoMuon);
1977  h_Criteria_NMuonTrack_den->SetBinContent(1 + 2 * i, 1 + 2 * j, CountRecoMuon);
1978  }
1979  }
1980  }
1981 
1982  h_Efficiency_ThetaXZ->Divide(h_ThetaXZ_num, h_ThetaXZ_den, 1, 1, "");
1983  h_Efficiency_ThetaYZ->Divide(h_ThetaYZ_num, h_ThetaYZ_den, 1, 1, "");
1985 
1989 
1992 
1997 
2002 
2007  }
2008  //========================================================================
2010 
2011 }
def analyze(root, level, gtrees, gbranches, doprint)
Definition: rootstat.py:69
intermediate_table::iterator iterator
void beginRun(const art::Run &run) override
unsigned int NumberTrajectoryPoints() const
Definition: MCParticle.h:218
const TLorentzVector & Position(const int i=0) const
Definition: MCParticle.h:219
int PdgCode() const
Definition: MCParticle.h:212
std::vector< sim::TrackIDE > HitToTrackIDEs(detinfo::DetectorClocksData const &clockData, recob::Hit const &hit) const
void FuncDistanceAndAngleBetweenTruthAndRecoTrack(const simb::MCParticle *&MCparticle, art::Ptr< recob::Track > Track, double &TempDistanceBetweenTruthAndRecoTrack, double &TempAngleBeetweenTruthAndRecoTrack)
std::string string
Definition: nybbler.cc:12
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
std::vector< TrackID > TrackIDs
const simb::MCParticle * TrackIdToParticle_P(int id) const
int Mother() const
Definition: MCParticle.h:213
constexpr T pow(T x)
Definition: pow.h:72
TH2D * h_MuonTrackStitching_FailedCriteria_CriteriaTwoTracks_Angle
Geometry information for a single TPC.
Definition: TPCGeo.h:38
struct vector vector
TH2D * h_MuonTrackStitching_MatchedCriteria_TrackResLeading_TrackResSecond
TH2D * h_Efficiency_ThetaXZ_SinThetaYZ_LeadingPlusSecond
Vector_t VertexDirection() const
Definition: Track.h:132
TH2D * h_MuonTrackStitching_FailedCriteria_Distance_Angle
bool get(SelectorBase const &, Handle< PROD > &result) const
Definition: DataViewImpl.h:606
std::set< const gar::rec::Track * > TrackList
Definition: TrackCreator.h:14
STL namespace.
TH2D * h_MuonTrackStitching_MatchedCriteria_TrackRes_Completeness
intermediate_table::const_iterator const_iterator
void truthMatcher(detinfo::DetectorClocksData const &clockData, std::vector< art::Ptr< recob::Hit >> AllHits, std::vector< art::Ptr< recob::Hit >> track_hits, const simb::MCParticle *&MCparticle, double &Purity, double &Completeness, double &TotalRecoEnergy)
Particle class.
static QStrList * l
Definition: config.cpp:1044
Definition: Run.h:17
int TrackId() const
Definition: MCParticle.h:210
TH2D * h_MuonTrackStitching_MatchedCriteria_CriteriaTwoTracks_Angle
TH2D * h_MuonTrackStitching_FailedCriteria_CompletenessSecondMuon_Angle
TH2D * h_MuonTrackStitching_FailedCriteria_CompletenessLeading_CompletenessSecond
bool insideFV(double vertex[4])
geo::TPCID FindTPCAtPosition(double const worldLoc[3]) const
Returns the ID of the TPC at specified location.
bool isRealData() const
TH2D * h_MuonTrackStitching_FailedCriteria_TrackRes_Completeness
bool getByLabel(std::string const &label, std::string const &instance, Handle< PROD > &result) const
Definition: DataViewImpl.h:633
TH2D * h_MuonTrackStitching_TrackResLeading_TrackResSecond
double Length(size_t p=0) const
Access to various track properties.
Definition: Track.h:167
TH2D * h_MuonTrackStitching_FailedCriteriaAndLeadingPlusSecondBad_TrackResLeading_TrackResSecond
#define DEFINE_ART_MODULE(klass)
Definition: ModuleMacros.h:67
void beginJob()
Definition: Breakpoints.cc:14
TH2D * h_Efficiency_ThetaXZ_ThetaYZ_DifferenceLeadingAndLeadingPlusSecond
T get(std::string const &key) const
Definition: ParameterSet.h:271
int CountBadLeadingMuonTrackAndLeadingPlusSecondBadTrackTooShort
TH2D * h_MuonTrackStitching_FailedCriteriaAndLeadingPlusSecondBad_TrackRes_Completeness
Point_t const & Vertex() const
Definition: Track.h:124
void processEff(const art::Event &evt, bool &isFiducial)
p
Definition: test.py:223
void FuncDistanceAndAngleBetweenTracks(art::Ptr< recob::Track > Track1, art::Ptr< recob::Track > Track2, double &TempDistanceBetweenTracks, double &TempAngleBetweenTracks, double &TempCriteriaTwoTracks)
The data type to uniquely identify a TPC.
Definition: geo_types.h:386
Definition of data types for geometry description.
Detector simulation of raw signals on wires.
TH2D * h_MuonTrackStitching_FailedCriteria_TrackResSecondMuon_Angle
const sim::ParticleList & ParticleList() const
double Vx(const int i=0) const
Definition: MCParticle.h:221
int ID() const
Definition: Track.h:198
TH2D * h_MuonTrackStitching_FailedCriteriaButLeadingPlusSecondGood_TrackResLeading_TrackResSecond
int CountBadLeadingMuonTrackAndLeadingPlusSecondBadTrackTooLong
Contains all timing reference information for the detector.
TH2D * h_MuonTrackStitching_FailedCriteriaAndLeadingPlusSecondBad_Distance_Angle
Vector_t EndDirection() const
Definition: Track.h:133
const TLorentzVector & Momentum(const int i=0) const
Definition: MCParticle.h:220
TH2D * h_MuonTrackStitching_CompletenessSecondMuon_Angle
#define MF_LOG_DEBUG(id)
Provides recob::Track data product.
void analyze(const art::Event &evt) override
double Vz(const int i=0) const
Definition: MCParticle.h:223
Point_t const & End() const
Definition: Track.h:125
list x
Definition: train.py:276
TrackCollectionProxyElement< TrackCollProxy > Track
Proxy to an element of a proxy collection of recob::Track objects.
Definition: Track.h:1036
TH2D * h_MuonTrackStitching_FailedCriteriaButLeadingPlusSecondGood_TrackRes_Completeness
art::ServiceHandle< geo::Geometry const > geom
double truthLength(const simb::MCParticle *MCparticle)
TH2D * h_MuonTrackStitching_FailedCriteria_TrackResLeading_TrackResSecond
TCEvent evt
Definition: DataStructs.cxx:7
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
Definition: Ptr.h:297
int CountBadLeadingMuonTrackAndOnlyOneMuonTrackTrackTooShort
LArSoft geometry interface.
Definition: ChannelGeo.h:16
int CountBadLeadingMuonTrackAndLeadingPlusSecondBadCompleteness
void LocalToWorld(const double *tpc, double *world) const
Transform point from local TPC frame to world frame.
Definition: TPCGeo.h:563
double Vy(const int i=0) const
Definition: MCParticle.h:222
TH2D * h_MuonTrackStitching_MatchedCriteria_Distance_Angle
QTextStream & endl(QTextStream &s)
Event finding and building.
TH2D * h_MuonTrackStitching_FailedCriteriaButLeadingPlusSecondGood_Distance_Angle
vertex reconstruction