LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
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
17 #include "larcore/CoreUtils/ServiceUtil.h"
18 #include "larcore/Geometry/Geometry.h"
19 #include "larcorealg/Geometry/TPCGeo.h"
20 #include "larcoreobj/SimpleTypesAndConstants/geo_types.h"
21 #include "lardata/DetectorInfoServices/DetectorClocksService.h"
22 #include "lardataobj/RecoBase/Track.h"
23 #include "larsim/MCCheater/BackTrackerService.h"
24 #include "larsim/MCCheater/ParticleInventoryService.h"
25 #include "nusimdata/SimulationBase/MCParticle.h"
26 
27 // Framework includes
28 #include "art/Framework/Core/EDAnalyzer.h"
29 #include "art/Framework/Core/ModuleMacros.h"
30 #include "art/Framework/Principal/Event.h"
31 #include "art/Framework/Principal/Handle.h"
32 #include "art/Framework/Services/Registry/ServiceHandle.h"
33 #include "art_root_io/TFileService.h"
34 #include "canvas/Persistency/Common/FindManyP.h"
35 #include "fhiclcpp/ParameterSet.h"
36 #include "messagefacility/MessageLogger/MessageLogger.h"
37 
38 // ROOT includes
39 #include "TH1.h"
40 #include "TH2.h"
41 #include "TStyle.h"
42 #include "TVector3.h"
43 
44 #define MAX_TRACKS 1000
45 using namespace std;
46 
47 //========================================================================
48 
49 namespace DUNE {
50 
51  class MuonTrackingEff : public art::EDAnalyzer {
52  public:
53  explicit MuonTrackingEff(fhicl::ParameterSet const& pset);
54 
55  private:
56  void beginJob() override;
57  void endJob() override;
58  void beginRun(const art::Run& run) override;
59  void analyze(const art::Event& evt) override;
60 
61  void processEff(const art::Event& evt, bool& isFiducial);
62 
63  void truthMatcher(detinfo::DetectorClocksData const& clockData,
64  std::vector<recob::Hit> const& AllHits,
65  std::vector<art::Ptr<recob::Hit>> track_hits,
66  const simb::MCParticle*& MCparticle,
67  double& Purity,
68  double& Completeness,
69  double& TotalRecoEnergy);
70 
71  void FuncDistanceAndAngleBetweenTracks(art::Ptr<recob::Track> Track1,
72  art::Ptr<recob::Track> Track2,
73  double& TempDistanceBetweenTracks,
74  double& TempAngleBetweenTracks,
75  double& TempCriteriaTwoTracks);
76 
77  void FuncDistanceAndAngleBetweenTruthAndRecoTrack(const simb::MCParticle*& MCparticle,
78  art::Ptr<recob::Track> Track,
79  double& TempDistanceBetweenTruthAndRecoTrack,
80  double& TempAngleBeetweenTruthAndRecoTrack);
81 
82  double truthLength(const simb::MCParticle* MCparticle);
83 
84  bool insideFV(double vertex[4]);
85 
86  void doEfficiencies();
87 
88  // the parameters we'll read from the .fcl
89  std::string fMCTruthModuleLabel;
90  std::string fTrackModuleLabel;
91  int fMuonPDGCode;
92 
93  // int MC_isCC;
94  // int MC_incoming_PDG;
95  // double MC_incoming_P[4];
96  double MCTruthMuonVertex[4];
97  // double MCTruthMuonStartMomentum[4];
98 
99  int MCTruthMuonID;
100  double MCTruthMuonMomentum;
101 
102  // double MCTruthMuonMomentumXZ=0;
103  // double MCTruthMuonMomentumYZ=0;
104  double MCTruthMuonThetaXZ = 0;
105  double MCTruthMuonThetaYZ = 0;
106 
107  //Counts to calculate efficiency and failed criteria
108  int EventCounter = 0;
109 
110  int CountMCTruthMuon = 0;
111  int CountRecoMuon = 0;
112 
113  int CountGoodLeadingMuonTrack = 0;
114  int CountNoRecoTracks = 0;
115  int CountNoMuonTracks = 0;
116  int CountBadLeadingMuonTrack = 0;
117  int CountCompleteness = 0;
118  int CountPurity = 0;
119  int CountTrackLengthTooShort = 0;
120  int CountTrackLengthTooLong = 0;
121 
122  int CountBadLeadingMuonTrackAndOnlyOneMuonTrack = 0;
123  int CountBadLeadingMuonTrackButLeadingPlusSecondGood = 0;
124  int CountBadLeadingMuonTrackAndLeadingPlusSecondBad = 0;
125 
126  int CountBadLeadingMuonTrackAndLeadingPlusSecondBadCompleteness = 0;
127  int CountBadLeadingMuonTrackAndLeadingPlusSecondBadPurity = 0;
128  int CountBadLeadingMuonTrackAndLeadingPlusSecondBadTrackTooShort = 0;
129  int CountBadLeadingMuonTrackAndLeadingPlusSecondBadTrackTooLong = 0;
130 
131  int CountBadLeadingMuonTrackAndOnlyOneMuonTrackCompleteness = 0;
132  int CountBadLeadingMuonTrackAndOnlyOneMuonTrackPurity = 0;
133  int CountBadLeadingMuonTrackAndOnlyOneMuonTrackTrackTooShort = 0;
134  int CountBadLeadingMuonTrackAndOnlyOneMuonTrackTrackTooLong = 0;
135 
136  double Criteria;
137  // int NMuonTracksTooShort=0;
138 
139  int GoodEvents1MuonTrack = 0;
140  int GoodEvents2MuonTrack = 0;
141  int GoodEvents3MuonTrack = 0;
142  int GoodEvents4OrMoreMuonTrack = 0;
143 
144  int BadEvents0MuonTrack = 0;
145  int BadEvents1MuonTrack = 0;
146  int BadEvents2MuonTrack = 0;
147  int BadEvents3MuonTrack = 0;
148  int BadEvents4OrMoreMuonTrack = 0;
149 
150  //TH1Ds
151  //Single Criteria and total reco energy
152  TH1D* h_Purity;
153  TH1D* h_Completeness;
154  TH1D* h_TrackRes;
155  TH1D* h_TotalRecoEnergy;
156  TH1D* h_TruthLength;
157  TH1D* h_VertexRes;
158  TH1D* h_DirectionRes;
159 
160  //Efficiency ThetaXZ
161  TH1D* h_Efficiency_ThetaXZ;
162  TH1D* h_ThetaXZ_den;
163  TH1D* h_ThetaXZ_num;
164 
165  //Efficiency ThetaYZ
166  TH1D* h_Efficiency_ThetaYZ;
167  TH1D* h_ThetaYZ_den;
168  TH1D* h_ThetaYZ_num;
169 
170  //Efficiency SinThetaYZ
171  TH1D* h_Efficiency_SinThetaYZ;
172  TH1D* h_SinThetaYZ_den;
173  TH1D* h_SinThetaYZ_num;
174 
175  //TH2Ds
176  //Efficiency ThetaXZ vs. ThetaYZ
177  TH2D* h_Efficiency_ThetaXZ_ThetaYZ;
178  TH2D* h_ThetaXZ_ThetaYZ_den;
179  TH2D* h_ThetaXZ_ThetaYZ_num;
180  TH2D* h_FailedReconstruction_ThetaXZ_ThetaYZ;
181 
182  //Efficiency ThetaXZ vs. SinThetaYZ
183  TH2D* h_Efficiency_ThetaXZ_SinThetaYZ;
184  TH2D* h_ThetaXZ_SinThetaYZ_den;
185  TH2D* h_ThetaXZ_SinThetaYZ_num;
186  TH2D* h_FailedReconstruction_ThetaXZ_SinThetaYZ;
187 
188  //Efficiency ThetaXZ vs. ThetaYZ after summing up leading plus second
189  TH2D* h_Efficiency_ThetaXZ_ThetaYZ_LeadingPlusSecond;
190  TH2D* h_ThetaXZ_ThetaYZ_LeadingPlusSecondOk;
191  TH2D* h_ThetaXZ_ThetaYZ_LeadingPlusSecondOk_num;
192 
193  //Difference in efficiency before and after summing up leading plus second: ThetaXZ vs. ThetaYZ
194  TH2D* h_Efficiency_ThetaXZ_ThetaYZ_DifferenceLeadingAndLeadingPlusSecond;
195 
196  //Efficiency ThetaXZ vs. SinThetaYZ after summing up leading plus second
197  TH2D* h_Efficiency_ThetaXZ_SinThetaYZ_LeadingPlusSecond;
198  TH2D* h_ThetaXZ_SinThetaYZ_LeadingPlusSecondOk;
199  TH2D* h_ThetaXZ_SinThetaYZ_LeadingPlusSecondOk_num;
200 
201  //Failed Criteria
202  TH2D* h_NoRecoTrackAtAll_ThetaXZ_SinThetaYZ;
203  TH2D* h_NoMuonTrack_ThetaXZ_SinThetaYZ;
204  TH2D* h_TrackTooShort_ThetaXZ_SinThetaYZ;
205  TH2D* h_TrackTooLong_ThetaXZ_SinThetaYZ;
206  TH2D* h_Completeness_ThetaXZ_SinThetaYZ;
207  TH2D* h_Purity_ThetaXZ_SinThetaYZ;
208 
209  //Criteria vs. NRecoTrack
210  TH2D* h_Criteria_NRecoTrack;
211  TH2D* h_Criteria_NRecoTrack_den;
212  TH2D* h_Criteria_NRecoTrack_num;
213 
214  //Criteria vs. NMuonTrack
215  TH2D* h_Criteria_NMuonTrack;
216  TH2D* h_Criteria_NMuonTrack_den;
217  TH2D* h_Criteria_NMuonTrack_num;
218 
219  //NoMuonTrack: Max length of no muon track vs. PDG code of that track (MC truth)
220  TH2D* h_NoMuonTrack_MaxTrackLength_PDGCode;
221 
222  //Stitching variables: all events
223  TH2D* h_MuonTrackStitching_TrackRes_Completeness;
224  TH2D* h_MuonTrackStitching_TrackResLeading_TrackResSecond;
225  TH2D* h_MuonTrackStitching_Distance_Angle;
226  TH2D* h_MuonTrackStitching_TrackResSecondMuon_Angle;
227  TH2D* h_MuonTrackStitching_CompletenessSecondMuon_Angle;
228  TH2D* h_MuonTrackStitching_CriteriaTwoTracks_Angle;
229 
230  //Stitching variables: bad events
231  TH2D* h_MuonTrackStitching_FailedCriteria_TrackRes_Completeness;
232  TH2D* h_MuonTrackStitching_FailedCriteria_TrackResLeading_TrackResSecond;
233  TH2D* h_MuonTrackStitching_FailedCriteria_CompletenessLeading_CompletenessSecond;
234  TH2D* h_MuonTrackStitching_FailedCriteria_Distance_Angle;
235  TH2D* h_MuonTrackStitching_FailedCriteria_TrackResSecondMuon_Angle;
236  TH2D* h_MuonTrackStitching_FailedCriteria_CompletenessSecondMuon_Angle;
237  TH2D* h_MuonTrackStitching_FailedCriteria_CriteriaTwoTracks_Angle;
238 
239  //Stitching variables: good events
240  TH2D* h_MuonTrackStitching_MatchedCriteria_TrackRes_Completeness;
241  TH2D* h_MuonTrackStitching_MatchedCriteria_TrackResLeading_TrackResSecond;
242  TH2D* h_MuonTrackStitching_MatchedCriteria_Distance_Angle;
243  TH2D* h_MuonTrackStitching_MatchedCriteria_CriteriaTwoTracks_Angle;
244 
245  //Stitching variables: bad events but leading plus second ok
246  TH2D* h_MuonTrackStitching_FailedCriteriaButLeadingPlusSecondGood_TrackRes_Completeness;
247  TH2D*
248  h_MuonTrackStitching_FailedCriteriaButLeadingPlusSecondGood_TrackResLeading_TrackResSecond;
249  TH2D* h_MuonTrackStitching_FailedCriteriaButLeadingPlusSecondGood_Distance_Angle;
250 
251  //Stitching variables: bad events and leading plus second not ok
252  TH2D* h_MuonTrackStitching_FailedCriteriaAndLeadingPlusSecondBad_TrackRes_Completeness;
253  TH2D* h_MuonTrackStitching_FailedCriteriaAndLeadingPlusSecondBad_TrackResLeading_TrackResSecond;
254  TH2D* h_MuonTrackStitching_FailedCriteriaAndLeadingPlusSecondBad_Distance_Angle;
255 
256  float fFidVolCutX;
257  float fFidVolCutY;
258  float fFidVolCutZ;
259 
260  float fFidVolXmin;
261  float fFidVolXmax;
262  float fFidVolYmin;
263  float fFidVolYmax;
264  float fFidVolZmin;
265  float fFidVolZmax;
266 
267  art::ServiceHandle<geo::Geometry const> geom;
268 
269  //My histograms
270  int NThetaXZBins = 36;
271  int ThetaXZBinMin = 0;
272  int ThetaXZBinMax = 360;
273 
274  int NThetaYZBins = 18;
275  int ThetaYZBinMin = -90;
276  int ThetaYZBinMax = 90;
277 
278  int NSinThetaYZBins = 18;
279  int SinThetaYZBinMin = -1;
280  int SinThetaYZBinMax = 1;
281 
282  int NCriteriaBins = 13;
283  double CriteriaBinMin = -0.25;
284  double CriteriaBinMax = 6.25;
285 
286  int NRecoTracksBins = 19;
287  double RecoTracksBinMin = -0.25;
288  double RecoTracksBinMax = 9.25;
289 
290  }; // class MuonTrackingEff
291 
292  //========================================================================
293  MuonTrackingEff::MuonTrackingEff(fhicl::ParameterSet const& p) : EDAnalyzer(p)
294  {
295  fMCTruthModuleLabel = p.get<std::string>("MCTruthModuleLabel");
296  fTrackModuleLabel = p.get<std::string>("TrackModuleLabel");
297  fMuonPDGCode = p.get<int>("MuonPDGCode");
298  fFidVolCutX = p.get<float>("FidVolCutX");
299  fFidVolCutY = p.get<float>("FidVolCutY");
300  fFidVolCutZ = p.get<float>("FidVolCutZ");
301  }
302  //========================================================================
303  void MuonTrackingEff::beginJob()
304  {
305  std::cout << "job begin..." << std::endl;
306  // Get geometry.
307  auto const* geo = lar::providerFrom<geo::Geometry>();
308  // Define histogram boundaries (cm).
309  // For now only draw cryostat=0.
310  double minx = 1e9;
311  double maxx = -1e9;
312  double miny = 1e9;
313  double maxy = -1e9;
314  double minz = 1e9;
315  double maxz = -1e9;
316  for (auto const& tpc : geo->Iterate<geo::TPCGeo>(geo::CryostatID{0})) {
317  auto const world = tpc.GetCenter();
318  if (minx > world.X() - tpc.HalfWidth()) minx = world.X() - tpc.HalfWidth();
319  if (maxx < world.X() + tpc.HalfWidth()) maxx = world.X() + tpc.HalfWidth();
320  if (miny > world.Y() - tpc.HalfHeight()) miny = world.Y() - tpc.HalfHeight();
321  if (maxy < world.Y() + tpc.HalfHeight()) maxy = world.Y() + tpc.HalfHeight();
322  if (minz > world.Z() - tpc.Length() / 2.) minz = world.Z() - tpc.Length() / 2.;
323  if (maxz < world.Z() + tpc.Length() / 2.) maxz = world.Z() + tpc.Length() / 2.;
324  }
325 
326  fFidVolXmin = minx + fFidVolCutX;
327  fFidVolXmax = maxx - fFidVolCutX;
328  fFidVolYmin = miny + fFidVolCutY;
329  fFidVolYmax = maxy - fFidVolCutY;
330  fFidVolZmin = minz + fFidVolCutZ;
331  fFidVolZmax = maxz - fFidVolCutZ;
332 
333  std::cout << "Fiducial volume:"
334  << "\n"
335  << fFidVolXmin << "\t< x <\t" << fFidVolXmax << "\n"
336  << fFidVolYmin << "\t< y <\t" << fFidVolYmax << "\n"
337  << fFidVolZmin << "\t< z <\t" << fFidVolZmax << "\n";
338 
339  art::ServiceHandle<art::TFileService const> tfs;
340 
341  //TH1D's
342  gStyle->SetTitleOffset(1.3, "Y");
343 
344  //Single Criteria and total reco energy
345  h_Purity =
346  tfs->make<TH1D>("h_Purity", "All events: Purity vs. # events; Purity; # events", 60, 0, 1.2);
347 
348  h_Completeness = tfs->make<TH1D>(
349  "h_Completeness", "All events: Completeness vs # events; Completeness; # events", 60, 0, 1.2);
350  h_Completeness->SetLineColor(kBlue);
351 
352  h_TrackRes =
353  tfs->make<TH1D>("h_TrackRes",
354  "All events: L_{reco}/L_{truth} vs. # events; L_{reco}/L_{truth}; # events;",
355  75,
356  0,
357  1.5);
358  h_TrackRes->SetLineColor(kRed);
359 
360  h_TotalRecoEnergy = tfs->make<TH1D>("h_TotalRecoEnergy",
361  "All events: Total reco energy (sum of all hits in all "
362  "tracks) vs. # events; E_{reco., tot.} [MeV]; # events",
363  100,
364  0,
365  1000);
366 
367  h_TruthLength = tfs->make<TH1D>(
368  "h_TruthLength",
369  "All events: truth muon length vs. # events; truth muon length [cm]; # events",
370  100,
371  0,
372  300);
373 
374  h_VertexRes = tfs->make<TH1D>(
375  "h_VertexRes",
376  "All events: Vertex residuals vs. # events; #Delta vertex_{truth-teco} [cm]; # events",
377  300,
378  0,
379  300);
380 
381  h_DirectionRes = tfs->make<TH1D>(
382  "h_DirectionRes",
383  "All events: Angular residuals vs. # events; #Delta#theta_{truth-reco} [#circ]; # events",
384  180,
385  0,
386  180);
387 
388  //Efficiency ThetaXZ
389  h_Efficiency_ThetaXZ =
390  tfs->make<TH1D>("h_Efficiency_ThetaXZ",
391  "Muon reco efficiency vs. #theta_{XZ}; #theta_{XZ} [#circ]; Efficiency",
392  NThetaXZBins,
393  ThetaXZBinMin,
394  ThetaXZBinMax);
395  h_ThetaXZ_den =
396  tfs->make<TH1D>("h_ThetaXZ_den",
397  "# generated muons vs. #theta_{XZ}; #theta_{XZ} [#circ]; # generated muons",
398  NThetaXZBins,
399  ThetaXZBinMin,
400  ThetaXZBinMax);
401  h_ThetaXZ_num = tfs->make<TH1D>(
402  "h_ThetaXZ_num",
403  "# reconstructed muons vs. #theta_{XZ}; #theta_{XZ} [#circ]; # reconstructed muons",
404  NThetaXZBins,
405  ThetaXZBinMin,
406  ThetaXZBinMax);
407 
408  //Efficiency ThetaYZ
409  h_Efficiency_ThetaYZ = tfs->make<TH1D>(
410  "h_Efficiency_ThetaYZ",
411  "Muon reco efficiency vs. #theta_{YZ}; #theta_{YZ} [#circ]; Muon reco efficiency",
412  NThetaYZBins,
413  ThetaYZBinMin,
414  ThetaYZBinMax);
415  ;
416  h_ThetaYZ_den =
417  tfs->make<TH1D>("h_ThetaYZ_den",
418  "# generated muons vs. #theta_{YZ}; #theta_{YZ} [#circ]; # generated muons",
419  NThetaYZBins,
420  ThetaYZBinMin,
421  ThetaYZBinMax);
422  h_ThetaYZ_num = tfs->make<TH1D>(
423  "h_ThetaYZ_num",
424  "# reconstructed muons vs. #theta_{YZ}; #theta_{YZ} [#circ]; # reconstructed muons",
425  NThetaYZBins,
426  ThetaYZBinMin,
427  ThetaYZBinMax);
428 
429  //Efficiency SinThetaYZ
430  h_Efficiency_SinThetaYZ = tfs->make<TH1D>(
431  "h_Efficiency_SinThetaYZ",
432  "Muon reco efficiency vs. sin(#theta_{YZ}); sin(#theta_{YZ}); Muon reco efficiency",
433  NSinThetaYZBins,
434  SinThetaYZBinMin,
435  SinThetaYZBinMax);
436  ;
437  h_SinThetaYZ_den =
438  tfs->make<TH1D>("h_SinThetaYZ_den",
439  "# generated muons vs. sin(#theta_{YZ}); sin(#theta_{YZ}); # generated muons",
440  NSinThetaYZBins,
441  SinThetaYZBinMin,
442  SinThetaYZBinMax);
443  h_SinThetaYZ_num = tfs->make<TH1D>(
444  "h_SinThetaYZ_num",
445  "# reconstructed muons vs. sin(#theta_{YZ}); sin(#theta_{YZ}); # reconstructed muons",
446  NSinThetaYZBins,
447  SinThetaYZBinMin,
448  SinThetaYZBinMax);
449 
450  h_Purity->Sumw2();
451  h_Completeness->Sumw2();
452  h_TrackRes->Sumw2();
453  h_TotalRecoEnergy->Sumw2();
454  h_TruthLength->Sumw2();
455  h_VertexRes->Sumw2();
456  h_DirectionRes->Sumw2();
457 
458  h_Efficiency_SinThetaYZ->Sumw2();
459  h_SinThetaYZ_den->Sumw2();
460  h_SinThetaYZ_num->Sumw2();
461 
462  h_Efficiency_ThetaXZ->Sumw2();
463  h_ThetaXZ_den->Sumw2();
464  h_ThetaXZ_num->Sumw2();
465 
466  h_Efficiency_ThetaYZ->Sumw2();
467  h_ThetaYZ_den->Sumw2();
468  h_ThetaYZ_num->Sumw2();
469 
470  //TH2D's
471  //Efficiency and Failed Reconstruction ThetaXZ vs. ThetaYZ
472  h_Efficiency_ThetaXZ_ThetaYZ = tfs->make<TH2D>(
473  "h_Efficiency_ThetaXZ_ThetaYZ",
474  "Muon reco efficiency: #theta_{XZ} vs. #theta_{YZ}; #theta_{XZ} [#circ]; #theta_{YZ} [#circ]",
475  NThetaXZBins,
476  ThetaXZBinMin,
477  ThetaXZBinMax,
478  NThetaYZBins,
479  ThetaYZBinMin,
480  ThetaYZBinMax);
481  h_Efficiency_ThetaXZ_ThetaYZ->SetOption("colz");
482 
483  h_ThetaXZ_ThetaYZ_den = tfs->make<TH2D>(
484  "h_ThetaXZ_ThetaYZ_den",
485  "# generated muons: #theta_{XZ} vs. #theta_{YZ}; #theta_{XZ} [#circ]; #theta_{YZ} [#circ]",
486  NThetaXZBins,
487  ThetaXZBinMin,
488  ThetaXZBinMax,
489  NThetaYZBins,
490  ThetaYZBinMin,
491  ThetaYZBinMax);
492  h_ThetaXZ_ThetaYZ_den->SetOption("colz");
493 
494  h_ThetaXZ_ThetaYZ_num = tfs->make<TH2D>("h_ThetaXZ_ThetaYZ_num",
495  "# reconstructed muons: #theta_{XZ} vs. #theta_{YZ}; "
496  "#theta_{XZ} [#circ]; #theta_{YZ} [#circ]",
497  NThetaXZBins,
498  ThetaXZBinMin,
499  ThetaXZBinMax,
500  NThetaYZBins,
501  ThetaYZBinMin,
502  ThetaYZBinMax);
503  h_ThetaXZ_ThetaYZ_num->SetOption("colz");
504 
505  h_FailedReconstruction_ThetaXZ_ThetaYZ =
506  tfs->make<TH2D>("h_FailedReconstruction_ThetaXZ_ThetaYZ",
507  "# failed reconstructions: #theta_{XZ} vs. #theta_{YZ}; #theta_{XZ} [#circ]; "
508  "#theta_{YZ} [#circ]",
509  NThetaXZBins,
510  ThetaXZBinMin,
511  ThetaXZBinMax,
512  NThetaYZBins,
513  ThetaYZBinMin,
514  ThetaYZBinMax);
515  h_FailedReconstruction_ThetaXZ_ThetaYZ->SetOption("colz");
516 
517  //Efficiency and Failed Reconstruction ThetaXZ vs. SinThetaYZ
518  h_Efficiency_ThetaXZ_SinThetaYZ =
519  tfs->make<TH2D>("h_Efficiency_ThetaXZ_SinThetaYZ",
520  "Muon reco efficiency: #theta_{XZ} vs. sin(#theta_{YZ}); #theta_{XZ} "
521  "[#circ]; sin(#theta_{YZ})",
522  NThetaXZBins,
523  ThetaXZBinMin,
524  ThetaXZBinMax,
525  NSinThetaYZBins,
526  SinThetaYZBinMin,
527  SinThetaYZBinMax);
528  h_Efficiency_ThetaXZ_SinThetaYZ->SetOption("colz");
529 
530  h_ThetaXZ_SinThetaYZ_den = tfs->make<TH2D>(
531  "h_ThetaXZ_SinThetaYZ_den",
532  "# generated muons: #theta_{XZ} vs. sin(#theta_{YZ}); #theta_{XZ} [#circ]; sin(#theta_{YZ})",
533  NThetaXZBins,
534  ThetaXZBinMin,
535  ThetaXZBinMax,
536  NSinThetaYZBins,
537  SinThetaYZBinMin,
538  SinThetaYZBinMax);
539  h_ThetaXZ_SinThetaYZ_den->SetOption("colz");
540 
541  h_ThetaXZ_SinThetaYZ_num =
542  tfs->make<TH2D>("h_ThetaXZ_SinThetaYZ_num",
543  "# reconstructed muons: #theta_{XZ} vs. sin(#theta_{YZ}); #theta_{XZ} "
544  "[#circ]; sin(#theta_{YZ}); #theta_{XZ} [#circ]; sin(#theta_{YZ})",
545  NThetaXZBins,
546  ThetaXZBinMin,
547  ThetaXZBinMax,
548  NSinThetaYZBins,
549  SinThetaYZBinMin,
550  SinThetaYZBinMax);
551  h_ThetaXZ_SinThetaYZ_num->SetOption("colz");
552 
553  h_FailedReconstruction_ThetaXZ_SinThetaYZ =
554  tfs->make<TH2D>("h_FailedReconstruction_ThetaXZ_SinThetaYZ",
555  "# failed reconstructions: #theta_{XZ} vs. sin(#theta_{YZ}); #theta_{XZ} "
556  "[#circ]; sin(#theta_{YZ})",
557  NThetaXZBins,
558  ThetaXZBinMin,
559  ThetaXZBinMax,
560  NSinThetaYZBins,
561  SinThetaYZBinMin,
562  SinThetaYZBinMax);
563  h_FailedReconstruction_ThetaXZ_SinThetaYZ->SetOption("colz");
564 
565  //Efficiency ThetaXZ vs. ThetaYZ after summing up leading plus second
566  h_Efficiency_ThetaXZ_ThetaYZ_LeadingPlusSecond =
567  tfs->make<TH2D>("h_Efficiency_ThetaXZ_ThetaYZ_LeadingPlusSecond",
568  "Muon reco efficiency after stitching: #theta_{XZ} vs. #theta_{YZ}; "
569  "#theta_{XZ} [#circ]; #theta_{YZ} [#circ]",
570  NThetaXZBins,
571  ThetaXZBinMin,
572  ThetaXZBinMax,
573  NThetaYZBins,
574  ThetaYZBinMin,
575  ThetaYZBinMax);
576  h_Efficiency_ThetaXZ_ThetaYZ_LeadingPlusSecond->SetOption("colz");
577 
578  h_ThetaXZ_ThetaYZ_LeadingPlusSecondOk =
579  tfs->make<TH2D>("h_ThetaXZ_ThetaYZ_LeadingPlusSecondOk",
580  "# reconstructed muons after stitching (failed before stitching): "
581  "#theta_{XZ} vs #theta_{YZ}; #theta_{XZ} [#circ]; #theta_{YZ} [#circ]",
582  NThetaXZBins,
583  ThetaXZBinMin,
584  ThetaXZBinMax,
585  NThetaYZBins,
586  ThetaYZBinMin,
587  ThetaYZBinMax);
588  h_ThetaXZ_ThetaYZ_LeadingPlusSecondOk->SetOption("colz");
589 
590  h_ThetaXZ_ThetaYZ_LeadingPlusSecondOk_num =
591  tfs->make<TH2D>("h_ThetaXZ_ThetaYZ_LeadingPlusSecondOk_num",
592  "# reconstructed muons after stitching: #theta_{XZ} vs. #theta_{YZ}; "
593  "#theta_{XZ} [#circ]; #theta_{YZ} [#circ]",
594  NThetaXZBins,
595  ThetaXZBinMin,
596  ThetaXZBinMax,
597  NThetaYZBins,
598  ThetaYZBinMin,
599  ThetaYZBinMax);
600  h_ThetaXZ_ThetaYZ_LeadingPlusSecondOk_num->SetOption("colz");
601 
602  //Efficiency ThetaXZ vs. SinThetaYZ after summing up leading plus second
603  h_Efficiency_ThetaXZ_SinThetaYZ_LeadingPlusSecond =
604  tfs->make<TH2D>("h_Efficiency_ThetaXZ_SinThetaYZ_LeadingPlusSecond",
605  "Muon reco efficiency after stitching: #theta_{XZ} vs. sin(#theta_{YZ}); "
606  "#theta_{XZ} [#circ]; sin(#theta_{YZ})",
607  NThetaXZBins,
608  ThetaXZBinMin,
609  ThetaXZBinMax,
610  NSinThetaYZBins,
611  SinThetaYZBinMin,
612  SinThetaYZBinMax);
613  h_Efficiency_ThetaXZ_SinThetaYZ_LeadingPlusSecond->SetOption("colz");
614 
615  h_ThetaXZ_SinThetaYZ_LeadingPlusSecondOk =
616  tfs->make<TH2D>("h_ThetaXZ_SinThetaYZ_LeadingPlusSecondOk",
617  "# reconstructed muons after stitching (failed before stitching): "
618  "#theta_{XZ} vs. sin(#theta_{YZ}); #theta_{XZ} [#circ]; sin(#theta_{YZ})",
619  NThetaXZBins,
620  ThetaXZBinMin,
621  ThetaXZBinMax,
622  NSinThetaYZBins,
623  SinThetaYZBinMin,
624  SinThetaYZBinMax);
625  h_ThetaXZ_SinThetaYZ_LeadingPlusSecondOk->SetOption("colz");
626 
627  h_ThetaXZ_SinThetaYZ_LeadingPlusSecondOk_num =
628  tfs->make<TH2D>("h_ThetaXZ_SinThetaYZ_LeadingPlusSecondOk_num",
629  "# reconstructed muons after stitching: #theta_{XZ} vs. sin(#theta_{YZ}); "
630  "#theta_{XZ} [#circ]; sin(#theta_{YZ})",
631  NThetaXZBins,
632  ThetaXZBinMin,
633  ThetaXZBinMax,
634  NSinThetaYZBins,
635  SinThetaYZBinMin,
636  SinThetaYZBinMax);
637  h_ThetaXZ_SinThetaYZ_LeadingPlusSecondOk_num->SetOption("colz");
638 
639  //Difference in efficiency before and after summing up leading plus second: ThetaXZ vs. ThetaYZ
640  h_Efficiency_ThetaXZ_ThetaYZ_DifferenceLeadingAndLeadingPlusSecond =
641  tfs->make<TH2D>("h_Efficiency_ThetaXZ_ThetaYZ_DifferenceLeadingAndLeadingPlusSecond",
642  "Muon reco efficiency: difference before and after stitching: #theta_{XZ} "
643  "vs. #theta_{YZ}; #theta_{XZ} [#circ]; #theta_{YZ} [#circ]",
644  NThetaXZBins,
645  ThetaXZBinMin,
646  ThetaXZBinMax,
647  NThetaYZBins,
648  ThetaYZBinMin,
649  ThetaYZBinMax);
650  h_Efficiency_ThetaXZ_ThetaYZ_DifferenceLeadingAndLeadingPlusSecond->SetOption("colz");
651 
652  //Failed Criteria
653  h_NoRecoTrackAtAll_ThetaXZ_SinThetaYZ =
654  tfs->make<TH2D>("h_NoRecoTrackAtAll_ThetaXZ_SinThetaYZ",
655  "# events with no reco track at all: #theta_{XZ} vs. sin(#theta_{YZ}); "
656  "#theta_{XZ} [#circ]; sin(#theta_{YZ})",
657  NThetaXZBins,
658  ThetaXZBinMin,
659  ThetaXZBinMax,
660  NSinThetaYZBins,
661  SinThetaYZBinMin,
662  SinThetaYZBinMax);
663  h_NoRecoTrackAtAll_ThetaXZ_SinThetaYZ->SetOption("colz");
664 
665  h_NoMuonTrack_ThetaXZ_SinThetaYZ =
666  tfs->make<TH2D>("h_NoMuonTrack_ThetaXZ_SinThetaYZ",
667  "# events with no muon track: #theta_{XZ} vs. sin(#theta_{YZ}); #theta_{XZ} "
668  "[#circ]; sin(#theta_{YZ})",
669  NThetaXZBins,
670  ThetaXZBinMin,
671  ThetaXZBinMax,
672  NSinThetaYZBins,
673  SinThetaYZBinMin,
674  SinThetaYZBinMax);
675  h_NoMuonTrack_ThetaXZ_SinThetaYZ->SetOption("colz");
676 
677  h_TrackTooShort_ThetaXZ_SinThetaYZ =
678  tfs->make<TH2D>("h_TrackTooShort_ThetaXZ_SinThetaYZ",
679  "# events with L_{reco}/L_{truth} < 75%: #theta_{XZ} vs. sin(#theta_{YZ}); "
680  "#theta_{XZ} [#circ]; sin(#theta_{YZ})",
681  NThetaXZBins,
682  ThetaXZBinMin,
683  ThetaXZBinMax,
684  NSinThetaYZBins,
685  SinThetaYZBinMin,
686  SinThetaYZBinMax);
687  h_TrackTooShort_ThetaXZ_SinThetaYZ->SetOption("colz");
688 
689  h_TrackTooLong_ThetaXZ_SinThetaYZ =
690  tfs->make<TH2D>("h_TrackTooLong_ThetaXZ_SinThetaYZ",
691  "#events with L_{reco}/L_{truth} > 125%: #theta_{XZ} vs. sin(#theta_{YZ}); "
692  "#theta_{XZ} [#circ]; sin(#theta_{YZ})",
693  NThetaXZBins,
694  ThetaXZBinMin,
695  ThetaXZBinMax,
696  NSinThetaYZBins,
697  SinThetaYZBinMin,
698  SinThetaYZBinMax);
699  h_TrackTooLong_ThetaXZ_SinThetaYZ->SetOption("colz");
700 
701  h_Completeness_ThetaXZ_SinThetaYZ =
702  tfs->make<TH2D>("h_Completeness_ThetaXZ_SinThetaYZ",
703  "# events with Completeness < 50%: #theta_{XZ} vs. sin(#theta_{YZ}); "
704  "#theta_{XZ} [#circ]; sin(#theta_{YZ})",
705  NThetaXZBins,
706  ThetaXZBinMin,
707  ThetaXZBinMax,
708  NSinThetaYZBins,
709  SinThetaYZBinMin,
710  SinThetaYZBinMax);
711  h_Completeness_ThetaXZ_SinThetaYZ->SetOption("colz");
712 
713  h_Purity_ThetaXZ_SinThetaYZ =
714  tfs->make<TH2D>("h_Purity_ThetaXZ_SinThetaYZ",
715  "# events with Purity < 50%: #theta_{XZ} vs. sin(#theta_{YZ}); #theta_{XZ} "
716  "[#circ]; sin(#theta_{YZ})",
717  NThetaXZBins,
718  ThetaXZBinMin,
719  ThetaXZBinMax,
720  NSinThetaYZBins,
721  SinThetaYZBinMin,
722  SinThetaYZBinMax);
723  h_Purity_ThetaXZ_SinThetaYZ->SetOption("colz");
724 
725  //Criteria vs. NRecoTrack
726  h_Criteria_NRecoTrack =
727  tfs->make<TH2D>("h_Criteria_NRecoTrack",
728  "Ratio: criteria vs. # reco tracks; Criteria; # reco tracks",
729  NCriteriaBins,
730  CriteriaBinMin,
731  CriteriaBinMax,
732  NRecoTracksBins,
733  RecoTracksBinMin,
734  RecoTracksBinMax);
735  h_Criteria_NRecoTrack->SetOption("colz");
736 
737  h_Criteria_NRecoTrack_num =
738  tfs->make<TH2D>("h_Criteria_NRecoTrack_num",
739  "# events: criteria vs. # reco tracks; Criteria; # reco tracks",
740  NCriteriaBins,
741  CriteriaBinMin,
742  CriteriaBinMax,
743  NRecoTracksBins,
744  RecoTracksBinMin,
745  RecoTracksBinMax);
746  h_Criteria_NRecoTrack_num->SetOption("colz");
747 
748  h_Criteria_NRecoTrack_den = tfs->make<TH2D>("h_Criteria_NRecoTrack_den",
749  "Divider histogram; Criteria; # reco tracks",
750  NCriteriaBins,
751  CriteriaBinMin,
752  CriteriaBinMax,
753  NRecoTracksBins,
754  RecoTracksBinMin,
755  RecoTracksBinMax);
756  h_Criteria_NRecoTrack_den->SetOption("colz");
757 
758  //Criteria vs. NMuonTrack
759  h_Criteria_NMuonTrack =
760  tfs->make<TH2D>("h_Criteria_NMuonTrack",
761  "Ratio: criteria vs. # muon tracks; Criteria; # muon tracks",
762  NCriteriaBins,
763  CriteriaBinMin,
764  CriteriaBinMax,
765  NRecoTracksBins,
766  RecoTracksBinMin,
767  RecoTracksBinMax);
768  h_Criteria_NMuonTrack->SetOption("colz");
769 
770  h_Criteria_NMuonTrack_num =
771  tfs->make<TH2D>("h_Criteria_NMuonTrack_num",
772  "# events: criteria vs. # muon tracks; Criteria; # muon tracks",
773  NCriteriaBins,
774  CriteriaBinMin,
775  CriteriaBinMax,
776  NRecoTracksBins,
777  RecoTracksBinMin,
778  RecoTracksBinMax);
779  h_Criteria_NMuonTrack_num->SetOption("colz");
780 
781  h_Criteria_NMuonTrack_den = tfs->make<TH2D>("h_Criteria_NMuonTrack_den",
782  "Divider histogram; Criteria; # muon tracks",
783  NCriteriaBins,
784  CriteriaBinMin,
785  CriteriaBinMax,
786  NRecoTracksBins,
787  RecoTracksBinMin,
788  RecoTracksBinMax);
789  h_Criteria_NMuonTrack_den->SetOption("colz");
790 
791  //NoMuonTrack: Max length of no muon track vs. PDG code of that track (MC truth)
792  h_NoMuonTrack_MaxTrackLength_PDGCode = tfs->make<TH2D>(
793  "h_NoMuonTrack_MaxTrackLength_PDGCode",
794  "Events with no muon track: L_{reco, max} vs. PDG Code; L_{reco} [cm]; PDG Code",
795  100,
796  0.,
797  200.,
798  200,
799  -100.,
800  100.);
801  h_NoMuonTrack_MaxTrackLength_PDGCode->SetOption("colz");
802 
803  //Stitching variables: all events
804  h_MuonTrackStitching_TrackRes_Completeness =
805  tfs->make<TH2D>("h_MuonTrackStitching_TrackRes_Completeness",
806  "All events: L_{reco}/L_{truth} (leading) vs. Completeness (leading); "
807  "L_{reco}/L_{truth} (leading); Completeness (leading)",
808  150,
809  0.,
810  1.5,
811  150,
812  0.,
813  1.5);
814  h_MuonTrackStitching_TrackRes_Completeness->SetOption("colz");
815 
816  h_MuonTrackStitching_TrackResLeading_TrackResSecond =
817  tfs->make<TH2D>("h_MuonTrackStitching_TrackResLeading_TrackResSecond",
818  "All events: L_{reco}/L_{truth} (leading) vs. L_{reco}/L_{truth} (second); "
819  "L_{reco}/L_{truth} (leading); L_{reco}/L_{truth} (second)",
820  150,
821  0.,
822  1.5,
823  150,
824  0.,
825  1.5);
826  h_MuonTrackStitching_TrackResLeading_TrackResSecond->SetOption("colz");
827 
828  h_MuonTrackStitching_Distance_Angle =
829  tfs->make<TH2D>("h_MuonTrackStitching_Distance_Angle",
830  "All events: distance vs. angle b/w leading and second muon track; Distance "
831  "[cm]; Angle [#circ]",
832  100,
833  0.,
834  100.,
835  100,
836  0.,
837  180.);
838  h_MuonTrackStitching_Distance_Angle->SetOption("colz");
839 
840  h_MuonTrackStitching_TrackResSecondMuon_Angle =
841  tfs->make<TH2D>("h_MuonTrackStitching_TrackResSecondMuon_Angle",
842  "All events: L_{reco}/L_{truth} (second) vs. angle; L_{reco}/L_{truth} "
843  "(second); Angle [#circ]",
844  150,
845  0.,
846  1.5,
847  180,
848  0,
849  180.);
850  h_MuonTrackStitching_TrackResSecondMuon_Angle->SetOption("colz");
851 
852  h_MuonTrackStitching_CompletenessSecondMuon_Angle = tfs->make<TH2D>(
853  "h_MuonTrackStitching_CompletenessSecondMuon_Angle",
854  "All events: Completeness (second) vs. angle; Completeness (second); Angle [#circ]",
855  120,
856  0.,
857  1.2,
858  180,
859  0,
860  180.);
861  h_MuonTrackStitching_CompletenessSecondMuon_Angle->SetOption("colz");
862 
863  h_MuonTrackStitching_CriteriaTwoTracks_Angle =
864  tfs->make<TH2D>("h_MuonTrackStitching_CriteriaTwoTracks_Angle",
865  "All events: CriteriaTwoTracks vs. angle; Criteria; Angle [#circ]",
866  7,
867  0.75,
868  4.25,
869  180,
870  0,
871  180.);
872  h_MuonTrackStitching_CriteriaTwoTracks_Angle->SetOption("colz");
873 
874  //Stitching variables: bad events
875  h_MuonTrackStitching_FailedCriteria_TrackRes_Completeness =
876  tfs->make<TH2D>("h_MuonTrackStitching_FailedCriteria_TrackRes_Completeness",
877  "Bad events: L_{reco}/L_{truth} (leading) vs. Completeness (leading); "
878  "L_{reco}/L_{truth} (leading); Completeness (leading)",
879  150,
880  0.,
881  1.5,
882  150,
883  0.,
884  1.5);
885  h_MuonTrackStitching_FailedCriteria_TrackRes_Completeness->SetOption("colz");
886 
887  h_MuonTrackStitching_FailedCriteria_Distance_Angle =
888  tfs->make<TH2D>("h_MuonTrackStitching_FailedCriteria_Distance_Angle",
889  "Bad events: distance vs. angle b/w leading and second muon track; Distance "
890  "[cm]; Angle [#circ]",
891  100,
892  0.,
893  100.,
894  100,
895  0.,
896  180.);
897  h_MuonTrackStitching_FailedCriteria_Distance_Angle->SetOption("colz");
898 
899  h_MuonTrackStitching_FailedCriteria_TrackResLeading_TrackResSecond =
900  tfs->make<TH2D>("h_MuonTrackStitching_FailedCriteria_TrackResLeading_TrackResSecond",
901  "Bad events: L_{reco}/L_{truth} (leading) vs. L_{reco}/L_{truth} (second); "
902  "L_{reco}/L_{truth} (leading); L_{reco}/L_{truth} (second)",
903  150,
904  0.,
905  1.5,
906  150,
907  0.,
908  1.5);
909  h_MuonTrackStitching_FailedCriteria_TrackResLeading_TrackResSecond->SetOption("colz");
910 
911  h_MuonTrackStitching_FailedCriteria_CompletenessLeading_CompletenessSecond =
912  tfs->make<TH2D>("*h_MuonTrackStitching_FailedCriteria_CompletenessLeading_CompletenessSecond",
913  "Bad events: Completeness (leading) vs. Completeness (second); Completeness "
914  "(leading); Completeness (second)",
915  150,
916  0.,
917  1.5,
918  150,
919  0.,
920  1.5);
921  h_MuonTrackStitching_FailedCriteria_CompletenessLeading_CompletenessSecond->SetOption("colz");
922 
923  h_MuonTrackStitching_FailedCriteria_TrackResSecondMuon_Angle =
924  tfs->make<TH2D>("h_MuonTrackStitching_FailedCriteria_TrackResSecondMuon_Angle",
925  "Bad events: L_{reco}/L_{truth} (second) vs. angle; L_{reco}/L_{truth} "
926  "(second); Angle [#circ]",
927  150,
928  0.,
929  1.5,
930  180,
931  0,
932  180.);
933  h_MuonTrackStitching_FailedCriteria_TrackResSecondMuon_Angle->SetOption("colz");
934 
935  h_MuonTrackStitching_FailedCriteria_CompletenessSecondMuon_Angle = tfs->make<TH2D>(
936  "h_MuonTrackStitching_FailedCriteria_CompletenessSecondMuon_Angle",
937  "Bad events: Completeness (second) vs. angle; Completeness (second); Angle [#circ]",
938  120,
939  0.,
940  1.2,
941  180,
942  0,
943  180.);
944  h_MuonTrackStitching_FailedCriteria_CompletenessSecondMuon_Angle->SetOption("colz");
945 
946  h_MuonTrackStitching_FailedCriteria_CriteriaTwoTracks_Angle =
947  tfs->make<TH2D>("h_MuonTrackStitching_FailedCriteria_CriteriaTwoTracks_Angle",
948  "Bad events: CriteriaTwoTracks vs. angle; Criteria; Angle [#circ]",
949  7,
950  0.75,
951  4.25,
952  180,
953  0,
954  180.);
955  h_MuonTrackStitching_FailedCriteria_CriteriaTwoTracks_Angle->SetOption("colz");
956 
957  //Stitching variables: good events
958  h_MuonTrackStitching_MatchedCriteria_TrackRes_Completeness =
959  tfs->make<TH2D>("h_MuonTrackStitching_MatchedCriteria_TrackRes_Completeness",
960  "Good events: L_{reco}/L_{truth} (leading) vs. Completeness (leading); "
961  "L_{reco}/L_{truth} (leading); Completeness (leading)",
962  150,
963  0.,
964  1.5,
965  150,
966  0.,
967  1.5);
968  h_MuonTrackStitching_MatchedCriteria_TrackRes_Completeness->SetOption("colz");
969 
970  h_MuonTrackStitching_MatchedCriteria_Distance_Angle =
971  tfs->make<TH2D>("h_MuonTrackStitching_MatchedCriteria_Distance_Angle",
972  "Good events: distance vs. angle b/w leading and second muon track; Distance "
973  "[cm]; Angle [#circ]",
974  100,
975  0.,
976  100.,
977  100,
978  0.,
979  180.);
980  h_MuonTrackStitching_MatchedCriteria_Distance_Angle->SetOption("colz");
981 
982  h_MuonTrackStitching_MatchedCriteria_TrackResLeading_TrackResSecond =
983  tfs->make<TH2D>("h_MuonTrackStitching_MatchedCriteria_TrackResLeading_TrackResSecond",
984  "Good events: L_{reco}/L_{truth} (leading) vs. L_{reco}/L_{truth} (second); "
985  "L_{reco}/L_{truth} (leading); L_{reco}/L_{truth} (second)",
986  150,
987  0.,
988  1.5,
989  150,
990  0.,
991  1.5);
992  h_MuonTrackStitching_MatchedCriteria_TrackResLeading_TrackResSecond->SetOption("colz");
993 
994  h_MuonTrackStitching_MatchedCriteria_CriteriaTwoTracks_Angle =
995  tfs->make<TH2D>("h_MuonTrackStitching_MatchedCriteria_CriteriaTwoTracks_Angle",
996  "Good events: CriteriaTwoTracks vs. angle b/w leading and second muon track; "
997  "Criteria; Angle [#circ]",
998  7,
999  0.75,
1000  4.25,
1001  100,
1002  0.,
1003  180.);
1004  h_MuonTrackStitching_MatchedCriteria_CriteriaTwoTracks_Angle->SetOption("colz");
1005 
1006  //Stitching variables: bad events but leading plus second ok
1007  h_MuonTrackStitching_FailedCriteriaButLeadingPlusSecondGood_TrackRes_Completeness =
1008  tfs->make<TH2D>(
1009  "h_MuonTrackStitching_FailedCriteriaButLeadingPlusSecondGood_TrackRes_Completeness",
1010  "Bad events but leading + second is good: L_{reco}/L_{truth} (leading) vs. Completeness "
1011  "(leading); L_{reco}/L_{truth} (leading); Completeness (leading)",
1012  150,
1013  0.,
1014  1.5,
1015  150,
1016  0.,
1017  1.5);
1018  h_MuonTrackStitching_FailedCriteriaButLeadingPlusSecondGood_TrackRes_Completeness->SetOption(
1019  "colz");
1020 
1021  h_MuonTrackStitching_FailedCriteriaButLeadingPlusSecondGood_Distance_Angle =
1022  tfs->make<TH2D>("h_MuonTrackStitching_FailedCriteriaButLeadingPlusSecondGood_Distance_Angle",
1023  "Bad events but leading + second is good: distance vs. angle b/w leading and "
1024  "second muon track; Distance [cm]; Angle [#circ]",
1025  100,
1026  0.,
1027  100.,
1028  100,
1029  0.,
1030  180.);
1031  h_MuonTrackStitching_FailedCriteriaButLeadingPlusSecondGood_Distance_Angle->SetOption("colz");
1032 
1033  h_MuonTrackStitching_FailedCriteriaButLeadingPlusSecondGood_TrackResLeading_TrackResSecond =
1034  tfs->make<TH2D>(
1035  "h_MuonTrackStitching_FailedCriteriaButLeadingPlusSecondGood_TrackResLeading_"
1036  "TrackResSecond",
1037  "Bad events but leading + second is good: L_{reco}/L_{truth} (leading) vs. "
1038  "L_{reco}/L_{truth} (second); L_{reco}/L_{truth} (leading); L_{reco}/L_{truth} (second)",
1039  150,
1040  0.,
1041  1.5,
1042  150,
1043  0.,
1044  1.5);
1045  h_MuonTrackStitching_FailedCriteriaButLeadingPlusSecondGood_TrackResLeading_TrackResSecond
1046  ->SetOption("colz");
1047 
1048  //Stitching variables: bad events and leading plus second not ok
1049  h_MuonTrackStitching_FailedCriteriaAndLeadingPlusSecondBad_TrackRes_Completeness =
1050  tfs->make<TH2D>(
1051  "h_MuonTrackStitching_FailedCriteriaAndLeadingPlusSecondBad_TrackRes_Completeness",
1052  "Bad events and leading + second is bad: L_{reco}/L_{truth} (leading) vs. Completeness "
1053  "(leading); L_{reco}/L_{truth} (leading); Completeness (leading)",
1054  150,
1055  0.,
1056  1.5,
1057  150,
1058  0.,
1059  1.5);
1060  h_MuonTrackStitching_FailedCriteriaAndLeadingPlusSecondBad_TrackRes_Completeness->SetOption(
1061  "colz");
1062 
1063  h_MuonTrackStitching_FailedCriteriaAndLeadingPlusSecondBad_Distance_Angle =
1064  tfs->make<TH2D>("h_MuonTrackStitching_FailedCriteriaAndLeadingPlusSecondBad_Distance_Angle",
1065  "Bad events and leading + second is bad: distance vs. angle b/w leading and "
1066  "second muon track; Distance [cm]; Angle [#circ]",
1067  100,
1068  0.,
1069  100.,
1070  100,
1071  0.,
1072  180.);
1073  h_MuonTrackStitching_FailedCriteriaAndLeadingPlusSecondBad_Distance_Angle->SetOption("colz");
1074 
1075  h_MuonTrackStitching_FailedCriteriaAndLeadingPlusSecondBad_TrackResLeading_TrackResSecond =
1076  tfs->make<TH2D>(
1077  "h_MuonTrackStitching_FailedCriteriaAndLeadingPlusSecondBad_TrackResLeading_TrackResSecond",
1078  "Bad events and leading + second is bad: L_{reco}/L_{truth} (leading) vs. "
1079  "L_{reco}/L_{truth} (second); L_{reco}/L_{truth} (leading); L_{reco}/L_{truth} (second)",
1080  150,
1081  0.,
1082  1.5,
1083  150,
1084  0.,
1085  1.5);
1086  h_MuonTrackStitching_FailedCriteriaAndLeadingPlusSecondBad_TrackResLeading_TrackResSecond
1087  ->SetOption("colz");
1088  }
1089  //========================================================================
1090  void MuonTrackingEff::endJob()
1091  {
1092  doEfficiencies();
1093  }
1094  //========================================================================
1095  void MuonTrackingEff::beginRun(const art::Run& /*run*/)
1096  {
1097  mf::LogInfo("MuonTrackingEff") << "begin run..." << std::endl;
1098  }
1099  //========================================================================
1100  void MuonTrackingEff::analyze(const art::Event& event)
1101  {
1102  if (event.isRealData()) return;
1103 
1104  bool isFiducial = false;
1105  processEff(event, isFiducial);
1106  }
1107  //========================================================================
1108  void MuonTrackingEff::processEff(const art::Event& event, bool& isFiducial)
1109  {
1110 
1111  EventCounter++;
1112  simb::MCParticle* MCTruthMuonParticle = NULL;
1113 
1114  art::ServiceHandle<cheat::ParticleInventoryService const> pi_serv;
1115  const sim::ParticleList& plist = pi_serv->ParticleList();
1116  simb::MCParticle* particle = 0;
1117 
1118  for (sim::ParticleList::const_iterator ipar = plist.begin(); ipar != plist.end(); ++ipar) {
1119  particle = ipar->second;
1120 
1121  if (particle->Mother() ==
1122  0) { //0=particle has no mother particle, 1=particle has a mother particle
1123  const TLorentzVector& positionStart = particle->Position(0);
1124  positionStart.GetXYZT(
1125  MCTruthMuonVertex); //MCTruthMuonVertex[0-2]=truth vertex, MCTruthMuonVertex[3]=t=0
1126  }
1127 
1128  if (particle->PdgCode() == fMuonPDGCode) { // Particle cannon muon
1129  MCTruthMuonParticle = particle;
1130  MCTruthMuonID = particle->TrackId();
1131  MCTruthMuonMomentum =
1132  sqrt(pow(particle->Momentum().Px(), 2) + pow(particle->Momentum().Py(), 2) +
1133  pow(particle->Momentum().Pz(), 2));
1134 
1135  if (particle->Momentum().Pz() >= 0 && particle->Momentum().Px() >= 0) {
1136  MCTruthMuonThetaXZ =
1137  (180.0 / 3.14159) * atan(particle->Momentum().Px() / particle->Momentum().Pz());
1138  }
1139  else if (particle->Momentum().Pz() < 0 && particle->Momentum().Px() >= 0) {
1140  MCTruthMuonThetaXZ =
1141  180.0 + (180.0 / 3.14159) * atan(particle->Momentum().Px() / particle->Momentum().Pz());
1142  }
1143  else if (particle->Momentum().Pz() < 0 && particle->Momentum().Px() < 0) {
1144  MCTruthMuonThetaXZ =
1145  180.0 + (180.0 / 3.14159) * atan(particle->Momentum().Px() / particle->Momentum().Pz());
1146  }
1147  else if (particle->Momentum().Pz() >= 0 && particle->Momentum().Px() < 0) {
1148  MCTruthMuonThetaXZ =
1149  360.0 + (180.0 / 3.14159) * atan(particle->Momentum().Px() / particle->Momentum().Pz());
1150  }
1151 
1152  MCTruthMuonThetaYZ =
1153  (180.0 / 3.14159) * asin(particle->Momentum().Py() / MCTruthMuonMomentum);
1154  }
1155  }
1156  double MCTruthLengthMuon = truthLength(MCTruthMuonParticle);
1157  h_TruthLength->Fill(MCTruthLengthMuon);
1158 
1159  //===================================================================
1160  //Saving denominator histograms
1161  //===================================================================
1162  isFiducial = insideFV(MCTruthMuonVertex);
1163  if (!isFiducial) return;
1164 
1165  //save events for Nucleon decay and particle cannon
1166  if (MCTruthMuonParticle) {
1167  h_ThetaXZ_den->Fill(MCTruthMuonThetaXZ);
1168  h_ThetaYZ_den->Fill(MCTruthMuonThetaYZ);
1169  h_SinThetaYZ_den->Fill(sin((3.14159 / 180.) * MCTruthMuonThetaYZ));
1170  h_ThetaXZ_ThetaYZ_den->Fill(MCTruthMuonThetaXZ, MCTruthMuonThetaYZ);
1171  h_ThetaXZ_SinThetaYZ_den->Fill(MCTruthMuonThetaXZ,
1172  sin((3.14159 / 180.) * MCTruthMuonThetaYZ));
1173  CountMCTruthMuon++;
1174  }
1175 
1176  //========================================================================
1177  // Reco stuff, once we have selected a MC Particle let's find out if there is a track associated
1178  //========================================================================
1179 
1180  int NMuonTracks = 0;
1181 
1182  art::Handle<std::vector<recob::Track>> TrackListHandle;
1183  if (!event.getByLabel(fTrackModuleLabel, TrackListHandle)) return;
1184  std::vector<art::Ptr<recob::Track>> TrackList;
1185  art::fill_ptr_vector(TrackList, TrackListHandle);
1186  int NRecoTracks = TrackList.size();
1187  art::FindManyP<recob::Hit> track_hits(TrackListHandle, event, fTrackModuleLabel);
1188  if (NRecoTracks == 0) {
1189  MF_LOG_DEBUG("MuonTrackingEff") << "There are no reco tracks... bye";
1190  std::cout << "There are no reco tracks! MCTruthMuonThetaXZ: " << std::endl;
1191 
1192  Criteria = 0.;
1193  h_FailedReconstruction_ThetaXZ_ThetaYZ->Fill(MCTruthMuonThetaXZ, MCTruthMuonThetaYZ);
1194  h_FailedReconstruction_ThetaXZ_SinThetaYZ->Fill(MCTruthMuonThetaXZ,
1195  sin((3.14159 / 180.) * MCTruthMuonThetaYZ));
1196  h_NoRecoTrackAtAll_ThetaXZ_SinThetaYZ->Fill(MCTruthMuonThetaXZ,
1197  sin((3.14159 / 180.) * MCTruthMuonThetaYZ));
1198  h_Criteria_NRecoTrack_num->Fill(Criteria, static_cast<double>(NRecoTracks));
1199  h_Criteria_NMuonTrack_num->Fill(Criteria, static_cast<double>(NMuonTracks));
1200  CountNoRecoTracks++;
1201  return;
1202  }
1203  MF_LOG_DEBUG("MuonTrackingEff") << "Found this many reco tracks " << NRecoTracks;
1204 
1205  //std::cout << "NRecoTracks: " << NRecoTracks << std::endl;
1206 
1207  double PurityLeadingMuon = 0.;
1208  double CompletenessLeadingMuon = 0.;
1209  double RecoLengthLeadingMuon = 0.;
1210  art::Ptr<recob::Track> TrackLeadingMuon;
1211 
1212  double RecoLengthSecondMuon = 0.;
1213  double CompletenessSecondMuon = 0.;
1214  double PuritySecondMuon = 0.;
1215  art::Ptr<recob::Track> TrackSecondMuon;
1216 
1217  double tmpTotalRecoEnergy = 0.;
1218 
1219  double MaxLengthNoRecoMuon = 0;
1220  int PDGCodeMaxLengthNoRecoMuon = 0;
1221 
1222  const simb::MCParticle* RecoMuonParticle = NULL;
1223 
1224  std::vector<art::Ptr<recob::Hit>> const& tmp_TrackHits = track_hits.at(0);
1225  std::vector<recob::Hit> const& AllHits = tmp_TrackHits[0].parentAs<std::vector>();
1226 
1227  auto const clockData =
1228  art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(event);
1229 
1230  // Loop over reco tracks
1231  for (int i = 0; i < NRecoTracks; i++) {
1232  art::Ptr<recob::Track> track = TrackList[i];
1233  std::vector<art::Ptr<recob::Hit>> TrackHits = track_hits.at(i);
1234  double tmpPurity = 0.;
1235  double tmpCompleteness = 0.;
1236  const simb::MCParticle* particle;
1237 
1238  truthMatcher(
1239  clockData, AllHits, TrackHits, particle, tmpPurity, tmpCompleteness, tmpTotalRecoEnergy);
1240 
1241  if (!particle) {
1242  std::cout << "ERROR: Truth matcher didn't find a particle!" << std::endl;
1243  continue;
1244  }
1245 
1246  if (track->Length() > MaxLengthNoRecoMuon && particle->PdgCode() != fMuonPDGCode &&
1247  particle->TrackId() != MCTruthMuonID) {
1248  MaxLengthNoRecoMuon = track->Length();
1249  PDGCodeMaxLengthNoRecoMuon = particle->PdgCode();
1250  }
1251  //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
1252  if ((particle->PdgCode() == fMuonPDGCode) && (particle->TrackId() == MCTruthMuonID) &&
1253  tmpCompleteness > 0 && tmpPurity > 0) {
1254 
1255  NMuonTracks++;
1256 
1257  if (NMuonTracks == 1) {
1258  CompletenessLeadingMuon = tmpCompleteness;
1259  PurityLeadingMuon = tmpPurity;
1260  RecoLengthLeadingMuon = track->Length();
1261  TrackLeadingMuon = track;
1262 
1263  RecoMuonParticle = particle;
1264  }
1265 
1266  if (NMuonTracks >= 2 && tmpCompleteness > CompletenessLeadingMuon) {
1267 
1268  CompletenessSecondMuon = CompletenessLeadingMuon;
1269  PuritySecondMuon = PurityLeadingMuon;
1270  RecoLengthSecondMuon = RecoLengthLeadingMuon;
1271  TrackSecondMuon = TrackLeadingMuon;
1272 
1273  CompletenessLeadingMuon = tmpCompleteness;
1274  PurityLeadingMuon = tmpPurity;
1275  RecoLengthLeadingMuon = track->Length();
1276  TrackLeadingMuon = track;
1277 
1278  RecoMuonParticle = particle;
1279  }
1280 
1281  else if (NMuonTracks >= 2 && tmpCompleteness < CompletenessLeadingMuon &&
1282  tmpCompleteness > CompletenessSecondMuon) {
1283  CompletenessSecondMuon = tmpCompleteness;
1284  PuritySecondMuon = tmpPurity;
1285  RecoLengthSecondMuon = track->Length();
1286  TrackSecondMuon = track;
1287  }
1288  }
1289  }
1290 
1291  h_TotalRecoEnergy->Fill(tmpTotalRecoEnergy);
1292 
1293  //Muon events
1294  //=========================================================
1295 
1296  if (RecoMuonParticle && MCTruthMuonParticle) {
1297  CountRecoMuon++;
1298  h_Purity->Fill(PurityLeadingMuon);
1299  h_Completeness->Fill(CompletenessLeadingMuon);
1300  h_TrackRes->Fill(RecoLengthLeadingMuon / MCTruthLengthMuon);
1301 
1302  std::cout << "TrackLeadingMuon->Vertex().X(): " << TrackLeadingMuon->Vertex().X()
1303  << std::endl;
1304  std::cout << "MCTruthMuonParticle->Vz(): " << MCTruthMuonParticle->Vz() << std::endl;
1305  std::cout << " " << std::endl;
1306 
1307  double DistanceBetweenTruthAndRecoTrack;
1308  double AngleBetweenTruthAndRecoTrack;
1309  FuncDistanceAndAngleBetweenTruthAndRecoTrack(RecoMuonParticle,
1310  TrackLeadingMuon,
1311  DistanceBetweenTruthAndRecoTrack,
1312  AngleBetweenTruthAndRecoTrack);
1313 
1314  h_VertexRes->Fill(DistanceBetweenTruthAndRecoTrack);
1315  h_DirectionRes->Fill(AngleBetweenTruthAndRecoTrack);
1316 
1317  h_MuonTrackStitching_TrackRes_Completeness->Fill(RecoLengthLeadingMuon / MCTruthLengthMuon,
1318  CompletenessLeadingMuon);
1319  if (NMuonTracks >= 2) {
1320  double DistanceBetweenTracks;
1321  double AngleBetweenTracks;
1322  double CriteriaTwoTracks;
1323 
1324  FuncDistanceAndAngleBetweenTracks(TrackSecondMuon,
1325  TrackLeadingMuon,
1326  DistanceBetweenTracks,
1327  AngleBetweenTracks,
1328  CriteriaTwoTracks);
1329 
1330  h_MuonTrackStitching_TrackResLeading_TrackResSecond->Fill(
1331  RecoLengthLeadingMuon / MCTruthLengthMuon, RecoLengthSecondMuon / MCTruthLengthMuon);
1332  h_MuonTrackStitching_Distance_Angle->Fill(DistanceBetweenTracks, AngleBetweenTracks);
1333  h_MuonTrackStitching_TrackResSecondMuon_Angle->Fill(
1334  RecoLengthSecondMuon / MCTruthLengthMuon, AngleBetweenTracks);
1335  h_MuonTrackStitching_CompletenessSecondMuon_Angle->Fill(CompletenessSecondMuon,
1336  AngleBetweenTracks);
1337  h_MuonTrackStitching_CriteriaTwoTracks_Angle->Fill(CriteriaTwoTracks, AngleBetweenTracks);
1338  }
1339 
1340  //Completeness
1341  if (CompletenessLeadingMuon < 0.5) {
1342  CountCompleteness++;
1343  Criteria = 2.;
1344  h_Completeness_ThetaXZ_SinThetaYZ->Fill(MCTruthMuonThetaXZ,
1345  sin((3.14159 / 180.) * MCTruthMuonThetaYZ));
1346  h_Criteria_NRecoTrack_num->Fill(Criteria, static_cast<double>(NRecoTracks));
1347  h_Criteria_NMuonTrack_num->Fill(Criteria, static_cast<double>(NMuonTracks));
1348  }
1349  //Purity
1350  if (PurityLeadingMuon < 0.5) {
1351  CountPurity++;
1352  Criteria = 3.;
1353  h_Purity_ThetaXZ_SinThetaYZ->Fill(MCTruthMuonThetaXZ,
1354  sin((3.14159 / 180.) * MCTruthMuonThetaYZ));
1355  h_Criteria_NRecoTrack_num->Fill(Criteria, static_cast<double>(NRecoTracks));
1356  h_Criteria_NMuonTrack_num->Fill(Criteria, static_cast<double>(NMuonTracks));
1357  }
1358  //Track too short
1359  if (RecoLengthLeadingMuon / MCTruthLengthMuon < 0.75) {
1360  CountTrackLengthTooShort++;
1361  Criteria = 4.;
1362  h_TrackTooShort_ThetaXZ_SinThetaYZ->Fill(MCTruthMuonThetaXZ,
1363  sin((3.14159 / 180.) * MCTruthMuonThetaYZ));
1364  h_Criteria_NRecoTrack_num->Fill(Criteria, static_cast<double>(NRecoTracks));
1365  h_Criteria_NMuonTrack_num->Fill(Criteria, static_cast<double>(NMuonTracks));
1366  }
1367  //Track too long
1368  if (RecoLengthLeadingMuon / MCTruthLengthMuon > 1.25) {
1369  CountTrackLengthTooLong++;
1370  Criteria = 5.;
1371  h_TrackTooLong_ThetaXZ_SinThetaYZ->Fill(MCTruthMuonThetaXZ,
1372  sin((3.14159 / 180.) * MCTruthMuonThetaYZ));
1373  h_Criteria_NRecoTrack_num->Fill(Criteria, static_cast<double>(NRecoTracks));
1374  h_Criteria_NMuonTrack_num->Fill(Criteria, static_cast<double>(NMuonTracks));
1375 
1376  /*std::cout << "Track too long!" << std::endl;
1377  art::ServiceHandle<cheat::BackTracker const> bt2;
1378  const sim::ParticleList& plist2 = bt2->ParticleList();
1379  simb::MCParticle *particle2=0;
1380 
1381  std::cout << "EventCounter: " << EventCounter << std::endl;
1382  for( sim::ParticleList::const_iterator ipar2 = plist2.begin(); ipar2!=plist2.end(); ++ipar2)
1383  {
1384  particle2 = ipar2->second;
1385  std::cout << "particle2->TrackId(): " << particle2->TrackId() << std::endl;
1386  std::cout << "particle2->PdgCode(): " << particle2->PdgCode() << std::endl;
1387  std::cout << "particle2->Momentum().P(): " << particle2->Momentum().P() << std::endl;
1388  std::cout << "tuthLength(particle2): " << truthLength(particle2) << std::endl;
1389  std::cout << "particle2->Position(): (x,y,z): " << particle2->Vx() << "\t" << particle2->Vy() << "\t" << particle2->Vz() << std::endl;
1390  std::cout << "particle2->Momentum(): (Px,Py,Pz): " << particle2->Momentum().Px() << "\t" << particle2->Momentum().Py() << "\t" << particle2->Momentum().Pz() << std::endl;
1391  std::cout << "particle2->Position().T(): " << particle2->Position().T() << std::endl;
1392  std::cout << "" << std::endl;
1393  }*/
1394  }
1395  //Reco failed at least one of the above criteria
1396  if (CompletenessLeadingMuon < 0.5 || PurityLeadingMuon < 0.5 ||
1397  RecoLengthLeadingMuon / MCTruthLengthMuon < 0.75 ||
1398  RecoLengthLeadingMuon / MCTruthLengthMuon > 1.25) {
1399  CountBadLeadingMuonTrack++;
1400  h_FailedReconstruction_ThetaXZ_ThetaYZ->Fill(MCTruthMuonThetaXZ, MCTruthMuonThetaYZ);
1401  h_FailedReconstruction_ThetaXZ_SinThetaYZ->Fill(MCTruthMuonThetaXZ,
1402  sin((3.14159 / 180.) * MCTruthMuonThetaYZ));
1403  h_MuonTrackStitching_FailedCriteria_TrackRes_Completeness->Fill(
1404  RecoLengthLeadingMuon / MCTruthLengthMuon, CompletenessLeadingMuon);
1405 
1406  if (NMuonTracks == 1) { BadEvents1MuonTrack++; }
1407  if (NMuonTracks == 2) { BadEvents2MuonTrack++; }
1408  if (NMuonTracks == 3) { BadEvents3MuonTrack++; }
1409  if (NMuonTracks == 4) { BadEvents4OrMoreMuonTrack++; }
1410 
1411  if (NMuonTracks >= 2) {
1412  double AngleBetweenTracks;
1413  double DistanceBetweenTracks;
1414  double CriteriaTwoTracks;
1415  FuncDistanceAndAngleBetweenTracks(TrackSecondMuon,
1416  TrackLeadingMuon,
1417  DistanceBetweenTracks,
1418  AngleBetweenTracks,
1419  CriteriaTwoTracks);
1420 
1421  if (AngleBetweenTracks > 160.) {
1422  std::cout << "EventCounter: " << EventCounter << std::endl;
1423  std::cout << "Angle b/w tracks: " << AngleBetweenTracks << std::endl;
1424  std::cout << "CriteriaTwoTracks: " << CriteriaTwoTracks << std::endl;
1425  std::cout << "CompletenessLeadingMuon: " << CompletenessLeadingMuon << std::endl;
1426  std::cout << "CompletenessSecondMuon: " << CompletenessSecondMuon << std::endl;
1427  std::cout << "PurityLeadingMuon: " << PurityLeadingMuon << std::endl;
1428  std::cout << "PuritySecondMuon: " << PuritySecondMuon << std::endl;
1429  std::cout << "TrackLeadingMuon: " << RecoLengthLeadingMuon / MCTruthLengthMuon
1430  << std::endl;
1431  std::cout << "TrackResSecondMuon: " << RecoLengthSecondMuon / MCTruthLengthMuon
1432  << std::endl;
1433  std::cout << "TrackID leading: " << TrackLeadingMuon->ID() << std::endl;
1434  std::cout << "TrackID second: " << TrackSecondMuon->ID() << std::endl;
1435  }
1436 
1437  h_MuonTrackStitching_FailedCriteria_Distance_Angle->Fill(DistanceBetweenTracks,
1438  AngleBetweenTracks);
1439  h_MuonTrackStitching_FailedCriteria_TrackResLeading_TrackResSecond->Fill(
1440  RecoLengthLeadingMuon / MCTruthLengthMuon, RecoLengthSecondMuon / MCTruthLengthMuon);
1441  h_MuonTrackStitching_FailedCriteria_CompletenessLeading_CompletenessSecond->Fill(
1442  CompletenessLeadingMuon, CompletenessSecondMuon);
1443  h_MuonTrackStitching_FailedCriteria_TrackResSecondMuon_Angle->Fill(
1444  RecoLengthSecondMuon / MCTruthLengthMuon, AngleBetweenTracks);
1445  h_MuonTrackStitching_FailedCriteria_CompletenessSecondMuon_Angle->Fill(
1446  CompletenessSecondMuon, AngleBetweenTracks);
1447  h_MuonTrackStitching_FailedCriteria_CriteriaTwoTracks_Angle->Fill(CriteriaTwoTracks,
1448  AngleBetweenTracks);
1449  if ((CompletenessLeadingMuon + CompletenessSecondMuon) >= 0.5 &&
1450  PurityLeadingMuon >= 0.5 && PuritySecondMuon >= 0.5 &&
1451  (RecoLengthLeadingMuon + RecoLengthSecondMuon) / MCTruthLengthMuon >= 0.75 &&
1452  (RecoLengthLeadingMuon + RecoLengthSecondMuon) / MCTruthLengthMuon <= 1.25) {
1453  CountBadLeadingMuonTrackButLeadingPlusSecondGood++;
1454  h_ThetaXZ_ThetaYZ_LeadingPlusSecondOk->Fill(MCTruthMuonThetaXZ, MCTruthMuonThetaYZ);
1455  h_ThetaXZ_SinThetaYZ_LeadingPlusSecondOk->Fill(
1456  MCTruthMuonThetaXZ, sin((3.14159 / 180.) * MCTruthMuonThetaYZ));
1457 
1458  h_MuonTrackStitching_FailedCriteriaButLeadingPlusSecondGood_TrackRes_Completeness->Fill(
1459  RecoLengthLeadingMuon / MCTruthLengthMuon, CompletenessLeadingMuon);
1460  h_MuonTrackStitching_FailedCriteriaButLeadingPlusSecondGood_TrackResLeading_TrackResSecond
1461  ->Fill(RecoLengthLeadingMuon / MCTruthLengthMuon,
1462  RecoLengthSecondMuon / MCTruthLengthMuon);
1463  h_MuonTrackStitching_FailedCriteriaButLeadingPlusSecondGood_Distance_Angle->Fill(
1464  DistanceBetweenTracks, AngleBetweenTracks);
1465  }
1466  if ((CompletenessLeadingMuon + CompletenessSecondMuon) < 0.5 || PurityLeadingMuon < 0.5 ||
1467  PuritySecondMuon < 0.5 ||
1468  (RecoLengthLeadingMuon + RecoLengthSecondMuon) / MCTruthLengthMuon < 0.75 ||
1469  (RecoLengthLeadingMuon + RecoLengthSecondMuon) / MCTruthLengthMuon > 1.25) {
1470  CountBadLeadingMuonTrackAndLeadingPlusSecondBad++;
1471 
1472  h_MuonTrackStitching_FailedCriteriaAndLeadingPlusSecondBad_TrackRes_Completeness->Fill(
1473  RecoLengthLeadingMuon / MCTruthLengthMuon, CompletenessLeadingMuon);
1474  h_MuonTrackStitching_FailedCriteriaAndLeadingPlusSecondBad_TrackResLeading_TrackResSecond
1475  ->Fill(RecoLengthLeadingMuon / MCTruthLengthMuon,
1476  RecoLengthSecondMuon / MCTruthLengthMuon);
1477  h_MuonTrackStitching_FailedCriteriaAndLeadingPlusSecondBad_Distance_Angle->Fill(
1478  DistanceBetweenTracks, AngleBetweenTracks);
1479  }
1480  if ((CompletenessLeadingMuon + CompletenessSecondMuon) < 0.5) {
1481  CountBadLeadingMuonTrackAndLeadingPlusSecondBadCompleteness++;
1482  }
1483  if (PurityLeadingMuon < 0.5 || PuritySecondMuon < 0.5) {
1484  CountBadLeadingMuonTrackAndLeadingPlusSecondBadPurity++;
1485  }
1486  if ((RecoLengthLeadingMuon + RecoLengthSecondMuon) / MCTruthLengthMuon < 0.75) {
1487  CountBadLeadingMuonTrackAndLeadingPlusSecondBadTrackTooShort++;
1488  }
1489  if ((RecoLengthLeadingMuon + RecoLengthSecondMuon) / MCTruthLengthMuon > 1.25) {
1490  CountBadLeadingMuonTrackAndLeadingPlusSecondBadTrackTooLong++;
1491  }
1492  }
1493  else if (NMuonTracks == 1) {
1494  CountBadLeadingMuonTrackAndOnlyOneMuonTrack++;
1495  if (CompletenessLeadingMuon < 0.5) {
1496  CountBadLeadingMuonTrackAndOnlyOneMuonTrackCompleteness++;
1497  }
1498  if (PurityLeadingMuon < 0.5) { CountBadLeadingMuonTrackAndOnlyOneMuonTrackPurity++; }
1499  if (RecoLengthLeadingMuon / MCTruthLengthMuon < 0.75) {
1500  CountBadLeadingMuonTrackAndOnlyOneMuonTrackTrackTooShort++;
1501  }
1502  if (RecoLengthLeadingMuon / MCTruthLengthMuon > 1.25) {
1503  CountBadLeadingMuonTrackAndOnlyOneMuonTrackTrackTooLong++;
1504  }
1505  }
1506  }
1507  //Reco succesful according to the above criteria
1508  if (CompletenessLeadingMuon >= 0.5 && PurityLeadingMuon >= 0.5 &&
1509  RecoLengthLeadingMuon / MCTruthLengthMuon >= 0.75 &&
1510  RecoLengthLeadingMuon / MCTruthLengthMuon <= 1.25) {
1511  CountGoodLeadingMuonTrack++;
1512  Criteria = 6.;
1513  h_ThetaXZ_num->Fill(MCTruthMuonThetaXZ);
1514  h_ThetaYZ_num->Fill(MCTruthMuonThetaYZ);
1515  h_SinThetaYZ_num->Fill(sin((3.14159 / 180.) * MCTruthMuonThetaYZ));
1516  h_ThetaXZ_ThetaYZ_num->Fill(MCTruthMuonThetaXZ, MCTruthMuonThetaYZ);
1517  h_ThetaXZ_SinThetaYZ_num->Fill(MCTruthMuonThetaXZ,
1518  sin((3.14159 / 180.) * MCTruthMuonThetaYZ));
1519  h_Criteria_NRecoTrack_num->Fill(Criteria, static_cast<double>(NRecoTracks));
1520  h_Criteria_NMuonTrack_num->Fill(Criteria, static_cast<double>(NMuonTracks));
1521 
1522  h_MuonTrackStitching_MatchedCriteria_TrackRes_Completeness->Fill(
1523  RecoLengthLeadingMuon / MCTruthLengthMuon, CompletenessLeadingMuon);
1524 
1525  if (NMuonTracks == 1) { GoodEvents1MuonTrack++; }
1526  if (NMuonTracks == 2) { GoodEvents2MuonTrack++; }
1527  if (NMuonTracks == 3) { GoodEvents3MuonTrack++; }
1528  if (NMuonTracks == 4) { GoodEvents4OrMoreMuonTrack++; }
1529 
1530  if (NMuonTracks >= 2) {
1531  h_MuonTrackStitching_MatchedCriteria_TrackResLeading_TrackResSecond->Fill(
1532  RecoLengthLeadingMuon / MCTruthLengthMuon, RecoLengthSecondMuon / MCTruthLengthMuon);
1533 
1534  double AngleBetweenTracks;
1535  double DistanceBetweenTracks;
1536  double CriteriaTwoTracks;
1537  FuncDistanceAndAngleBetweenTracks(TrackSecondMuon,
1538  TrackLeadingMuon,
1539  DistanceBetweenTracks,
1540  AngleBetweenTracks,
1541  CriteriaTwoTracks);
1542  h_MuonTrackStitching_MatchedCriteria_Distance_Angle->Fill(DistanceBetweenTracks,
1543  AngleBetweenTracks);
1544  h_MuonTrackStitching_MatchedCriteria_CriteriaTwoTracks_Angle->Fill(CriteriaTwoTracks,
1545  AngleBetweenTracks);
1546  }
1547  }
1548  }
1549  //No muon track
1550  if (!RecoMuonParticle && MCTruthMuonParticle) {
1551  CountNoMuonTracks++;
1552  BadEvents0MuonTrack++;
1553  Criteria = 1.;
1554  h_NoMuonTrack_ThetaXZ_SinThetaYZ->Fill(MCTruthMuonThetaXZ,
1555  sin((3.14159 / 180.) * MCTruthMuonThetaYZ));
1556  h_FailedReconstruction_ThetaXZ_ThetaYZ->Fill(MCTruthMuonThetaXZ, MCTruthMuonThetaYZ);
1557  h_FailedReconstruction_ThetaXZ_SinThetaYZ->Fill(MCTruthMuonThetaXZ,
1558  sin((3.14159 / 180.) * MCTruthMuonThetaYZ));
1559  h_Criteria_NRecoTrack_num->Fill(Criteria, static_cast<double>(NRecoTracks));
1560  h_Criteria_NMuonTrack_num->Fill(Criteria, static_cast<double>(NMuonTracks));
1561  h_NoMuonTrack_MaxTrackLength_PDGCode->Fill(MaxLengthNoRecoMuon,
1562  static_cast<double>(PDGCodeMaxLengthNoRecoMuon));
1563  }
1564  }
1565  //========================================================================
1566  void MuonTrackingEff::truthMatcher(detinfo::DetectorClocksData const& clockData,
1567  std::vector<recob::Hit> const& AllHits,
1568  std::vector<art::Ptr<recob::Hit>> track_hits,
1569  const simb::MCParticle*& MCparticle,
1570  double& Purity,
1571  double& Completeness,
1572  double& TotalRecoEnergy)
1573  {
1574  art::ServiceHandle<cheat::BackTrackerService const> bt_serv;
1575  art::ServiceHandle<cheat::ParticleInventoryService const> pi_serv;
1576  std::map<int, double> trkID_E; // map that connects TrackID and energy for
1577  // each hit <trackID, energy>
1578  for (size_t j = 0; j < track_hits.size(); ++j) {
1579  art::Ptr<recob::Hit> hit = track_hits[j];
1580  std::vector<sim::TrackIDE> TrackIDs =
1581  bt_serv->HitToTrackIDEs(clockData,
1582  hit); // TrackIDE contains TrackID, energy and energyFrac. A hit can
1583  // have several TrackIDs (so this hit is associated with multiple
1584  // MC truth track IDs (EM shower IDs are negative). If a hit ahs
1585  // multiple trackIDs, "energyFrac" contains the fraction of the
1586  // energy of for each ID compared to the total energy of the hit.
1587  // "energy" contains only the energy associated with the specific
1588  // ID in that case. This requires MC truth info!
1589  for (size_t k = 0; k < TrackIDs.size(); k++) {
1590  trkID_E[TrackIDs[k].trackID] +=
1591  TrackIDs[k].energy; // sum up the energy for each TrackID and store
1592  // <TrackID, energy> in "TrkID_E"
1593  }
1594  }
1595 
1596  double max_E = -999.0;
1597  double TotalEnergyTrack = 0.0;
1598  int TrackID = -999;
1599  double PartialEnergyTrackID =
1600  0.0; // amount of energy deposited by the particle that deposited more
1601  // energy... tomato potato... blabla
1605  if (!trkID_E.size()) {
1606  MCparticle = 0;
1607  return; // Ghost track???
1608  }
1609  for (std::map<int, double>::iterator ii = trkID_E.begin(); ii != trkID_E.end();
1610  ++ii) { // trkID_E contains the trackID (first) and corresponding
1611  // energy (second) for a specific track, summed up over all
1612  // events. here looping over all trekID_E's
1613  TotalEnergyTrack += ii->second; // and summing up the energy of all hits
1614  // in the track (TotalEnergyTrack)
1615  if ((ii->second) > max_E) { // looking for the trakID with the highest energy in the
1616  // track. this is PartialEnergyTrackID and max_E then
1617  PartialEnergyTrackID = ii->second;
1618  max_E = ii->second;
1619  TrackID = ii->first; // saving trackID of the ID with the highest energy
1620  // contribution in the track to assign it to
1621  // MCparticle later
1622  }
1623  }
1624  MCparticle = pi_serv->TrackIdToParticle_P(TrackID);
1625  // In the current simulation, we do not save EM Shower daughters in GEANT.
1626  // But we do save the energy deposition in TrackIDEs. If the energy
1627  // deposition is from a particle that is the daughter of an EM particle, the
1628  // negative of the parent track ID is saved in TrackIDE for the daughter
1629  // particle we don't want to track gammas or any other EM activity
1630  if (TrackID < 0) { return; }
1631 
1632  Purity = PartialEnergyTrackID / TotalEnergyTrack;
1633 
1634  // completeness
1635  TotalRecoEnergy = 0;
1636  for (recob::Hit const& hit : AllHits) {
1637  std::vector<sim::TrackIDE> TrackIDs = bt_serv->HitToTrackIDEs(clockData, hit);
1638  for (size_t l = 0; l < TrackIDs.size(); ++l) { // and over all track IDs of the hits
1639  if (TrackIDs[l].trackID == TrackID)
1640  TotalRecoEnergy += TrackIDs[l].energy; // and sum up the energy fraction of all hits
1641  // that correspond ot the saved trackID
1642  }
1643  }
1644  Completeness = PartialEnergyTrackID / TotalRecoEnergy;
1645  }
1646 
1647  void MuonTrackingEff::FuncDistanceAndAngleBetweenTracks(art::Ptr<recob::Track> Track1,
1648  art::Ptr<recob::Track> Track2,
1649  double& TempDistanceBetweenTracks,
1650  double& TempAngleBetweenTracks,
1651  double& TempCriteriaTwoTracks)
1652  {
1653 
1654  TempDistanceBetweenTracks = sqrt(pow(Track1->End().X() - Track2->Vertex().X(), 2) +
1655  pow(Track1->End().Y() - Track2->Vertex().Y(), 2) +
1656  pow(Track1->End().Z() - Track2->Vertex().Z(), 2));
1657  TempAngleBetweenTracks = (180.0 / 3.14159) * Track1->EndDirection<TVector3>().Angle(
1658  Track2->VertexDirection<TVector3>());
1659  TempCriteriaTwoTracks = 1.;
1660 
1661  if (TempDistanceBetweenTracks > sqrt(pow(Track1->End().X() - Track2->End().X(), 2) +
1662  pow(Track1->End().Y() - Track2->End().Y(), 2) +
1663  pow(Track1->End().Z() - Track2->End().Z(), 2))) {
1664  TempDistanceBetweenTracks = sqrt(pow(Track1->End().X() - Track2->End().X(), 2) +
1665  pow(Track1->End().Y() - Track2->End().Y(), 2) +
1666  pow(Track1->End().Z() - Track2->End().Z(), 2));
1667  TempAngleBetweenTracks = 180. - (180.0 / 3.14159) * Track1->EndDirection<TVector3>().Angle(
1668  Track2->EndDirection<TVector3>());
1669  TempCriteriaTwoTracks = 2.;
1670  }
1671 
1672  if (TempDistanceBetweenTracks > sqrt(pow(Track1->Vertex().X() - Track2->End().X(), 2) +
1673  pow(Track1->Vertex().Y() - Track2->End().Y(), 2) +
1674  pow(Track1->Vertex().Z() - Track2->End().Z(), 2))) {
1675  TempDistanceBetweenTracks = sqrt(pow(Track1->Vertex().X() - Track2->End().X(), 2) +
1676  pow(Track1->Vertex().Y() - Track2->End().Y(), 2) +
1677  pow(Track1->Vertex().Z() - Track2->End().Z(), 2));
1678  TempAngleBetweenTracks = (180.0 / 3.14159) * Track1->VertexDirection<TVector3>().Angle(
1679  Track2->EndDirection<TVector3>());
1680  TempCriteriaTwoTracks = 3.;
1681  }
1682 
1683  if (TempDistanceBetweenTracks > sqrt(pow(Track1->Vertex().X() - Track2->Vertex().X(), 2) +
1684  pow(Track1->Vertex().Y() - Track2->Vertex().Y(), 2) +
1685  pow(Track1->Vertex().Z() - Track2->Vertex().Z(), 2))) {
1686  TempDistanceBetweenTracks = sqrt(pow(Track1->Vertex().X() - Track2->Vertex().X(), 2) +
1687  pow(Track1->Vertex().Y() - Track2->Vertex().Y(), 2) +
1688  pow(Track1->Vertex().Z() - Track2->Vertex().Z(), 2));
1689  TempAngleBetweenTracks = 180. - (180.0 / 3.14159) * Track1->VertexDirection<TVector3>().Angle(
1690  Track2->VertexDirection<TVector3>());
1691  TempCriteriaTwoTracks = 4.;
1692  }
1693  }
1694 
1695  void MuonTrackingEff::FuncDistanceAndAngleBetweenTruthAndRecoTrack(
1696  const simb::MCParticle*& MCparticle,
1697  art::Ptr<recob::Track> Track,
1698  double& TempDistanceBetweenTruthAndRecoTrack,
1699  double& TempAngleBeetweenTruthAndRecoTrack)
1700  {
1701  TempDistanceBetweenTruthAndRecoTrack = sqrt(pow(Track->Vertex().X() - MCparticle->Vx(), 2) +
1702  pow(Track->Vertex().Y() - MCparticle->Vy(), 2) +
1703  pow(Track->Vertex().Z() - MCparticle->Vz(), 2));
1704 
1705  TempAngleBeetweenTruthAndRecoTrack = 0;
1706  }
1707 
1708  //========================================================================
1709  double MuonTrackingEff::truthLength(const simb::MCParticle* MCparticle)
1710  {
1711  //calculate the truth length considering only the part that is inside the TPC
1712  //Base on a piece of code from dune/TrackingAna/TrackingEfficiency_module.cc
1713 
1714  if (!MCparticle) return -999.0;
1715  int numberTrajectoryPoints = MCparticle->NumberTrajectoryPoints();
1716  std::vector<double> TPCLengthHits(numberTrajectoryPoints, 0);
1717  int FirstHit = 0, LastHit = 0;
1718  double TPCLength = 0.0;
1719  bool BeenInVolume = false;
1720 
1721  for (unsigned int MCHit = 0; MCHit < TPCLengthHits.size(); ++MCHit) {
1722  if (MCHit != 0)
1723  TPCLengthHits[MCHit] = sqrt(pow((MCparticle->Vx(MCHit - 1) - MCparticle->Vx(MCHit)), 2) +
1724  pow((MCparticle->Vy(MCHit - 1) - MCparticle->Vy(MCHit)), 2) +
1725  pow((MCparticle->Vz(MCHit - 1) - MCparticle->Vz(MCHit)), 2));
1726  geo::TPCID tpcid =
1727  geom->FindTPCAtPosition(geo::vect::toPoint(MCparticle->Position(MCHit).Vect()));
1728  if (tpcid.isValid) {
1729  LastHit = MCHit;
1730  if (!BeenInVolume) {
1731  BeenInVolume = true;
1732  FirstHit = MCHit;
1733  }
1734  }
1735  }
1736  for (int Hit = FirstHit + 1; Hit <= LastHit; ++Hit)
1737  TPCLength += TPCLengthHits[Hit];
1738  return TPCLength;
1739  }
1740  //========================================================================
1741  bool MuonTrackingEff::insideFV(double vertex[4])
1742  {
1743  double x = vertex[0];
1744  double y = vertex[1];
1745  double z = vertex[2];
1746 
1747  return x > fFidVolXmin && x < fFidVolXmax && y > fFidVolYmin && y < fFidVolYmax &&
1748  z > fFidVolZmin && z < fFidVolZmax;
1749  }
1750  //========================================================================
1751  void MuonTrackingEff::doEfficiencies()
1752  {
1753  std::cout << std::endl;
1754 
1755  std::cout << "EventCounter: "
1756  << "\t" << EventCounter << std::endl;
1757 
1758  std::cout << "CountMCTruthMuon: "
1759  << "\t" << CountMCTruthMuon << " = "
1760  << 100 * static_cast<double>(CountMCTruthMuon) / static_cast<double>(EventCounter)
1761  << "%" << std::endl;
1762 
1763  std::cout << "CountGoodLeadingMuonTrack (=good events): "
1764  << "\t" << CountGoodLeadingMuonTrack << "/" << CountMCTruthMuon << " = "
1765  << 100 * static_cast<double>(CountGoodLeadingMuonTrack) /
1766  static_cast<double>(CountMCTruthMuon)
1767  << "%" << std::endl;
1768 
1769  std::cout << "CountBadLeadingMuonTrack+CountNoRecoTracks+CountNoMuonTracks (=bad events): "
1770  << "\t" << CountBadLeadingMuonTrack + CountNoRecoTracks + CountNoMuonTracks << " = "
1771  << 100 *
1772  static_cast<double>(CountBadLeadingMuonTrack + CountNoRecoTracks +
1773  CountNoMuonTracks) /
1774  static_cast<double>(CountMCTruthMuon)
1775  << "%" << std::endl;
1776 
1777  std::cout << "CountNoRecoTracks+CountNoMuonTracks: "
1778  << "\t" << CountNoRecoTracks + CountNoMuonTracks << " = "
1779  << 100 * static_cast<double>(CountNoRecoTracks + CountNoMuonTracks) /
1780  static_cast<double>(CountMCTruthMuon)
1781  << "%" << std::endl;
1782 
1783  std::cout << "CountTrackLengthTooShort: "
1784  << "\t" << CountTrackLengthTooShort << " = "
1785  << 100 * static_cast<double>(CountTrackLengthTooShort) /
1786  static_cast<double>(CountMCTruthMuon)
1787  << "%" << std::endl;
1788 
1789  std::cout << "CountCompleteness: "
1790  << "\t" << CountCompleteness << " = "
1791  << 100 * static_cast<double>(CountCompleteness) /
1792  static_cast<double>(CountMCTruthMuon)
1793  << "%" << std::endl;
1794 
1795  std::cout << "CountTrackLengthTooLong: "
1796  << "\t" << CountTrackLengthTooLong << " = "
1797  << 100 * static_cast<double>(CountTrackLengthTooLong) /
1798  static_cast<double>(CountMCTruthMuon)
1799  << "%" << std::endl;
1800 
1801  std::cout << "CountPurity: "
1802  << "\t" << CountPurity << " = "
1803  << 100 * static_cast<double>(CountPurity) / static_cast<double>(CountMCTruthMuon)
1804  << "%" << std::endl;
1805 
1806  std::cout << std::endl;
1807 
1808  std::cout << "GoodLeadingMuonTrack+CountBadLeadingMuonTrackButLeadingPlusSecondGood (=good "
1809  "events after stitching): "
1810  << "\t"
1811  << CountGoodLeadingMuonTrack + CountBadLeadingMuonTrackButLeadingPlusSecondGood << "/"
1812  << CountMCTruthMuon << " = "
1813  << 100 *
1814  static_cast<double>(CountGoodLeadingMuonTrack +
1815  CountBadLeadingMuonTrackButLeadingPlusSecondGood) /
1816  static_cast<double>(CountMCTruthMuon)
1817  << "%" << std::endl;
1818 
1819  std::cout << "CountBadLeadingMuonTrack+CountNoRecoTracks+CountNoMuonTracks-"
1820  "CountBadLeadingMuonTrackButLeadingPlusSecondGood (=bad events after stitching) : "
1821  << "\t"
1822  << CountBadLeadingMuonTrack + CountNoRecoTracks + CountNoMuonTracks -
1823  CountBadLeadingMuonTrackButLeadingPlusSecondGood
1824  << " = "
1825  << 100 *
1826  static_cast<double>(CountBadLeadingMuonTrack + CountNoRecoTracks +
1827  CountNoMuonTracks -
1828  CountBadLeadingMuonTrackButLeadingPlusSecondGood) /
1829  static_cast<double>(CountMCTruthMuon)
1830  << "%" << std::endl;
1831 
1832  std::cout << std::endl;
1833 
1834  std::cout << "CountBadLeadingMuonTrackButLeadingPlusSecondGood: "
1835  << "\t" << CountBadLeadingMuonTrackButLeadingPlusSecondGood << " = "
1836  << 100 * static_cast<double>(CountBadLeadingMuonTrackButLeadingPlusSecondGood) /
1837  static_cast<double>(CountMCTruthMuon)
1838  << "%" << std::endl;
1839 
1840  std::cout << "CountBadLeadingMuonTrackAndOnlyOneMuonTrack: "
1841  << "\t" << CountBadLeadingMuonTrackAndOnlyOneMuonTrack << " = "
1842  << 100 * static_cast<double>(CountBadLeadingMuonTrackAndOnlyOneMuonTrack) /
1843  static_cast<double>(CountMCTruthMuon)
1844  << "%" << std::endl;
1845 
1846  std::cout << "CountBadLeadingMuonTrackAndLeadingPlusSecondBad: "
1847  << "\t" << CountBadLeadingMuonTrackAndLeadingPlusSecondBad << " = "
1848  << 100 * static_cast<double>(CountBadLeadingMuonTrackAndLeadingPlusSecondBad) /
1849  static_cast<double>(CountMCTruthMuon)
1850  << "%" << std::endl;
1851 
1852  std::cout << "CountNoRecoTracks: "
1853  << "\t" << CountNoRecoTracks << " = "
1854  << 100 * static_cast<double>(CountNoRecoTracks) /
1855  static_cast<double>(CountMCTruthMuon)
1856  << "%" << std::endl;
1857 
1858  std::cout << "CountNoMuonTracks: "
1859  << "\t" << CountNoMuonTracks << " = "
1860  << 100 * static_cast<double>(CountNoMuonTracks) /
1861  static_cast<double>(CountMCTruthMuon)
1862  << "%" << std::endl;
1863 
1864  std::cout << std::endl;
1865 
1866  std::cout << "CountBadLeadingMuonTrackAndOnlyOneMuonTrackCompleteness: "
1867  << CountBadLeadingMuonTrackAndOnlyOneMuonTrackCompleteness << std::endl;
1868  std::cout << "CountBadLeadingMuonTrackAndOnlyOneMuonTrackPurity: "
1869  << CountBadLeadingMuonTrackAndOnlyOneMuonTrackPurity << std::endl;
1870  std::cout << "CountBadLeadingMuonTrackAndOnlyOneMuonTrackTrackTooShort: "
1871  << CountBadLeadingMuonTrackAndOnlyOneMuonTrackTrackTooShort << std::endl;
1872  std::cout << "CountBadLeadingMuonTrackAndOnlyOneMuonTrackTrackTooLong: "
1873  << CountBadLeadingMuonTrackAndOnlyOneMuonTrackTrackTooLong << std::endl;
1874 
1875  std::cout << std::endl;
1876 
1877  std::cout << "CountBadLeadingMuonTrackAndLeadingPlusSecondBadCompleteness: "
1878  << "\t" << CountBadLeadingMuonTrackAndLeadingPlusSecondBadCompleteness << " = "
1879  << 100 *
1880  static_cast<double>(
1881  CountBadLeadingMuonTrackAndLeadingPlusSecondBadCompleteness) /
1882  static_cast<double>(CountMCTruthMuon)
1883  << "%" << std::endl;
1884 
1885  std::cout << "CountBadLeadingMuonTrackAndLeadingPlusSecondBadPurity: "
1886  << "\t" << CountBadLeadingMuonTrackAndLeadingPlusSecondBadPurity << " = "
1887  << 100 * static_cast<double>(CountBadLeadingMuonTrackAndLeadingPlusSecondBadPurity) /
1888  static_cast<double>(CountMCTruthMuon)
1889  << "%" << std::endl;
1890 
1891  std::cout << "CountBadLeadingMuonTrackAndLeadingPlusSecondBadTrackTooShort: "
1892  << "\t" << CountBadLeadingMuonTrackAndLeadingPlusSecondBadTrackTooShort << " = "
1893  << 100 *
1894  static_cast<double>(
1895  CountBadLeadingMuonTrackAndLeadingPlusSecondBadTrackTooShort) /
1896  static_cast<double>(CountMCTruthMuon)
1897  << "%" << std::endl;
1898 
1899  std::cout << "CountBadLeadingMuonTrackAndLeadingPlusSecondBadTrackTooLong: "
1900  << "\t" << CountBadLeadingMuonTrackAndLeadingPlusSecondBadTrackTooLong << " = "
1901  << 100 *
1902  static_cast<double>(
1903  CountBadLeadingMuonTrackAndLeadingPlusSecondBadTrackTooLong) /
1904  static_cast<double>(CountMCTruthMuon)
1905  << "%" << std::endl;
1906 
1907  std::cout << std::endl;
1908 
1909  std::cout << "GoodEvents1MuonTrack: " << GoodEvents1MuonTrack << std::endl;
1910  std::cout << "GoodEvents2MuonTrack: " << GoodEvents2MuonTrack << std::endl;
1911  std::cout << "GoodEvents3MuonTrack: " << GoodEvents3MuonTrack << std::endl;
1912  std::cout << "GoodEvents4OrMoreMuonTrack: " << GoodEvents4OrMoreMuonTrack << std::endl;
1913 
1914  std::cout << "BadEvents0MuonTrack: " << BadEvents0MuonTrack << std::endl;
1915  std::cout << "BadEvents1MuonTrack: " << BadEvents1MuonTrack << std::endl;
1916  std::cout << "BadEvents2MuonTrack: " << BadEvents2MuonTrack << std::endl;
1917  std::cout << "BadEvents3MuonTrack: " << BadEvents3MuonTrack << std::endl;
1918  std::cout << "BadEvents4OrMoreMuonTrack: " << BadEvents4OrMoreMuonTrack << std::endl;
1919 
1920  art::ServiceHandle<art::TFileService const> tfs;
1921 
1922  for (int i = 0; i < (NCriteriaBins + 1) / 2; ++i) {
1923  for (int j = 0; j < (NRecoTracksBins + 1) / 2; ++j) {
1924  if (i == 0) {
1925  h_Criteria_NRecoTrack_den->SetBinContent(1 + 2 * i, 1 + 2 * j, CountNoRecoTracks);
1926  h_Criteria_NMuonTrack_den->SetBinContent(1 + 2 * i, 1 + 2 * j, CountNoRecoTracks);
1927  }
1928  else if (i == 1) {
1929  h_Criteria_NRecoTrack_den->SetBinContent(1 + 2 * i, 1 + 2 * j, CountNoMuonTracks);
1930  h_Criteria_NMuonTrack_den->SetBinContent(1 + 2 * i, 1 + 2 * j, CountNoMuonTracks);
1931  }
1932  else if (i == 2) {
1933  h_Criteria_NRecoTrack_den->SetBinContent(1 + 2 * i, 1 + 2 * j, CountCompleteness);
1934  h_Criteria_NMuonTrack_den->SetBinContent(1 + 2 * i, 1 + 2 * j, CountCompleteness);
1935  }
1936  else if (i == 3) {
1937  h_Criteria_NRecoTrack_den->SetBinContent(1 + 2 * i, 1 + 2 * j, CountPurity);
1938  h_Criteria_NMuonTrack_den->SetBinContent(1 + 2 * i, 1 + 2 * j, CountPurity);
1939  }
1940  else if (i == 4) {
1941  h_Criteria_NRecoTrack_den->SetBinContent(1 + 2 * i, 1 + 2 * j, CountTrackLengthTooShort);
1942  h_Criteria_NMuonTrack_den->SetBinContent(1 + 2 * i, 1 + 2 * j, CountTrackLengthTooShort);
1943  }
1944  else if (i == 5) {
1945  h_Criteria_NRecoTrack_den->SetBinContent(1 + 2 * i, 1 + 2 * j, CountTrackLengthTooLong);
1946  h_Criteria_NMuonTrack_den->SetBinContent(1 + 2 * i, 1 + 2 * j, CountTrackLengthTooLong);
1947  }
1948  else if (i == 6) {
1949  h_Criteria_NRecoTrack_den->SetBinContent(1 + 2 * i, 1 + 2 * j, CountRecoMuon);
1950  h_Criteria_NMuonTrack_den->SetBinContent(1 + 2 * i, 1 + 2 * j, CountRecoMuon);
1951  }
1952  }
1953  }
1954 
1955  h_Efficiency_ThetaXZ->Divide(h_ThetaXZ_num, h_ThetaXZ_den, 1, 1, "");
1956  h_Efficiency_ThetaYZ->Divide(h_ThetaYZ_num, h_ThetaYZ_den, 1, 1, "");
1957  h_Efficiency_SinThetaYZ->Divide(h_SinThetaYZ_num, h_SinThetaYZ_den, 1, 1, "");
1958 
1959  h_Efficiency_ThetaXZ_ThetaYZ->Divide(h_ThetaXZ_ThetaYZ_num, h_ThetaXZ_ThetaYZ_den, 1, 1, "");
1960  h_Efficiency_ThetaXZ_SinThetaYZ->Divide(
1961  h_ThetaXZ_SinThetaYZ_num, h_ThetaXZ_SinThetaYZ_den, 1, 1, "");
1962 
1963  h_Criteria_NRecoTrack->Divide(h_Criteria_NRecoTrack_num, h_Criteria_NRecoTrack_den, 1, 1, "");
1964  h_Criteria_NMuonTrack->Divide(h_Criteria_NMuonTrack_num, h_Criteria_NMuonTrack_den, 1, 1, "");
1965 
1966  h_ThetaXZ_ThetaYZ_LeadingPlusSecondOk_num->Add(h_ThetaXZ_ThetaYZ_num, 1);
1967  h_ThetaXZ_ThetaYZ_LeadingPlusSecondOk_num->Add(h_ThetaXZ_ThetaYZ_LeadingPlusSecondOk, 1);
1968  h_Efficiency_ThetaXZ_ThetaYZ_LeadingPlusSecond->Divide(
1969  h_ThetaXZ_ThetaYZ_LeadingPlusSecondOk_num, h_ThetaXZ_ThetaYZ_den, 1, 1, "");
1970 
1971  h_ThetaXZ_SinThetaYZ_LeadingPlusSecondOk_num->Add(h_ThetaXZ_SinThetaYZ_num, 1);
1972  h_ThetaXZ_SinThetaYZ_LeadingPlusSecondOk_num->Add(h_ThetaXZ_SinThetaYZ_LeadingPlusSecondOk, 1);
1973  h_Efficiency_ThetaXZ_SinThetaYZ_LeadingPlusSecond->Divide(
1974  h_ThetaXZ_SinThetaYZ_LeadingPlusSecondOk_num, h_ThetaXZ_SinThetaYZ_den, 1, 1, "");
1975 
1976  h_Efficiency_ThetaXZ_ThetaYZ_DifferenceLeadingAndLeadingPlusSecond->Add(
1977  h_Efficiency_ThetaXZ_ThetaYZ_LeadingPlusSecond, 1);
1978  h_Efficiency_ThetaXZ_ThetaYZ_DifferenceLeadingAndLeadingPlusSecond->Add(
1979  h_Efficiency_ThetaXZ_ThetaYZ, -1);
1980  }
1981  //========================================================================
1982  DEFINE_ART_MODULE(MuonTrackingEff)
1983 
1984 }
x
Float_t x
Definition: compare.C:6
iterator
intermediate_table::iterator iterator
Definition: intermediate_table.cc:27
DUNE::MuonTrackingEff::beginRun
void beginRun(const art::Run &run) override
Definition: MuonTrackingEff_module.cc:1095
simb::MCParticle::NumberTrajectoryPoints
unsigned int NumberTrajectoryPoints() const
Definition: MCParticle.h:219
simb::MCParticle::Position
const TLorentzVector & Position(const int i=0) const
Definition: MCParticle.h:220
simb::MCParticle::PdgCode
int PdgCode() const
Definition: MCParticle.h:213
ServiceUtil.h
Utilities related to art service access.
DUNE::MuonTrackingEff::truthMatcher
void truthMatcher(detinfo::DetectorClocksData const &clockData, std::vector< recob::Hit > const &AllHits, std::vector< art::Ptr< recob::Hit >> track_hits, const simb::MCParticle *&MCparticle, double &Purity, double &Completeness, double &TotalRecoEnergy)
Definition: MuonTrackingEff_module.cc:1566
cheat::BackTrackerService::HitToTrackIDEs
std::vector< sim::TrackIDE > HitToTrackIDEs(detinfo::DetectorClocksData const &clockData, recob::Hit const &hit) const
Definition: BackTrackerService.cc:163
DUNE::MuonTrackingEff::FuncDistanceAndAngleBetweenTruthAndRecoTrack
void FuncDistanceAndAngleBetweenTruthAndRecoTrack(const simb::MCParticle *&MCparticle, art::Ptr< recob::Track > Track, double &TempDistanceBetweenTruthAndRecoTrack, double &TempAngleBeetweenTruthAndRecoTrack)
Definition: MuonTrackingEff_module.cc:1695
mf::LogInfo
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
Definition: MessageLogger.h:209
TrackIDs
std::vector< TrackID > TrackIDs
Definition: ClusteringValidation_module.cc:60
cheat::ParticleInventoryService::TrackIdToParticle_P
const simb::MCParticle * TrackIdToParticle_P(int id) const
Definition: ParticleInventoryService.cc:154
simb::MCParticle::Mother
int Mother() const
Definition: MCParticle.h:214
y
Float_t y
Definition: compare.C:6
sim::ParticleList::const_iterator
list_type::const_iterator const_iterator
Definition: ParticleList.h:132
z
Double_t z
Definition: plot.C:276
DUNE::MuonTrackingEff::h_MuonTrackStitching_FailedCriteria_CriteriaTwoTracks_Angle
TH2D * h_MuonTrackStitching_FailedCriteria_CriteriaTwoTracks_Angle
Definition: MuonTrackingEff_module.cc:237
geo::TPCGeo
Geometry information for a single TPC.
Definition: TPCGeo.h:36
DUNE::MuonTrackingEff::h_MuonTrackStitching_MatchedCriteria_TrackResLeading_TrackResSecond
TH2D * h_MuonTrackStitching_MatchedCriteria_TrackResLeading_TrackResSecond
Definition: MuonTrackingEff_module.cc:241
geo::vect::toPoint
::geo::Point_t toPoint(Point const &p)
Convert the specified point into a geo::Point_t.
Definition: geo_vectors_utils.h:1584
DUNE::MuonTrackingEff::h_Efficiency_ThetaXZ_SinThetaYZ_LeadingPlusSecond
TH2D * h_Efficiency_ThetaXZ_SinThetaYZ_LeadingPlusSecond
Definition: MuonTrackingEff_module.cc:197
recob::Track::VertexDirection
Vector_t VertexDirection() const
Access to track direction at different points.
Definition: Track.h:166
DUNE::MuonTrackingEff::h_MuonTrackStitching_FailedCriteria_Distance_Angle
TH2D * h_MuonTrackStitching_FailedCriteria_Distance_Angle
Definition: MuonTrackingEff_module.cc:234
std
STL namespace.
DUNE::MuonTrackingEff::h_MuonTrackStitching_MatchedCriteria_TrackRes_Completeness
TH2D * h_MuonTrackStitching_MatchedCriteria_TrackRes_Completeness
Definition: MuonTrackingEff_module.cc:240
art::Event::isRealData
bool isRealData() const
Definition: Event.cc:53
MCParticle.h
Particle class.
art::Run
Definition: Run.h:37
simb::MCParticle::TrackId
int TrackId() const
Definition: MCParticle.h:211
DUNE::MuonTrackingEff::h_MuonTrackStitching_MatchedCriteria_CriteriaTwoTracks_Angle
TH2D * h_MuonTrackStitching_MatchedCriteria_CriteriaTwoTracks_Angle
Definition: MuonTrackingEff_module.cc:243
DUNE::MuonTrackingEff::h_MuonTrackStitching_FailedCriteria_CompletenessSecondMuon_Angle
TH2D * h_MuonTrackStitching_FailedCriteria_CompletenessSecondMuon_Angle
Definition: MuonTrackingEff_module.cc:236
DUNE::MuonTrackingEff::h_MuonTrackStitching_FailedCriteria_CompletenessLeading_CompletenessSecond
TH2D * h_MuonTrackStitching_FailedCriteria_CompletenessLeading_CompletenessSecond
Definition: MuonTrackingEff_module.cc:233
DUNE::MuonTrackingEff::insideFV
bool insideFV(double vertex[4])
Definition: MuonTrackingEff_module.cc:1741
lar::dump::vector
auto vector(Vector const &v)
Returns a manipulator which will print the specified array.
Definition: DumpUtils.h:289
DUNE::MuonTrackingEff::h_MuonTrackStitching_FailedCriteria_TrackRes_Completeness
TH2D * h_MuonTrackStitching_FailedCriteria_TrackRes_Completeness
Definition: MuonTrackingEff_module.cc:231
geo::GeometryCore::FindTPCAtPosition
TPCID FindTPCAtPosition(Point_t const &point) const
Returns the ID of the TPC at specified location.
Definition: GeometryCore.cxx:326
DUNE::MuonTrackingEff::h_MuonTrackStitching_TrackResLeading_TrackResSecond
TH2D * h_MuonTrackStitching_TrackResLeading_TrackResSecond
Definition: MuonTrackingEff_module.cc:224
recob::Track::Length
double Length(size_t p=0) const
Access to various track properties.
Definition: Track.h:207
DUNE::MuonTrackingEff::h_MuonTrackStitching_FailedCriteriaAndLeadingPlusSecondBad_TrackResLeading_TrackResSecond
TH2D * h_MuonTrackStitching_FailedCriteriaAndLeadingPlusSecondBad_TrackResLeading_TrackResSecond
Definition: MuonTrackingEff_module.cc:253
DEFINE_ART_MODULE
#define DEFINE_ART_MODULE(klass)
Definition: ModuleMacros.h:65
breakpoints::beginJob
void beginJob()
Definition: Breakpoints.cc:14
DUNE::MuonTrackingEff::h_Efficiency_ThetaXZ_ThetaYZ_DifferenceLeadingAndLeadingPlusSecond
TH2D * h_Efficiency_ThetaXZ_ThetaYZ_DifferenceLeadingAndLeadingPlusSecond
Definition: MuonTrackingEff_module.cc:194
fhicl::ParameterSet::get
T get(std::string const &key) const
Definition: ParameterSet.h:314
DUNE::MuonTrackingEff::CountBadLeadingMuonTrackAndLeadingPlusSecondBadTrackTooShort
int CountBadLeadingMuonTrackAndLeadingPlusSecondBadTrackTooShort
Definition: MuonTrackingEff_module.cc:128
Event.h
DUNE::MuonTrackingEff::h_MuonTrackStitching_FailedCriteriaAndLeadingPlusSecondBad_TrackRes_Completeness
TH2D * h_MuonTrackStitching_FailedCriteriaAndLeadingPlusSecondBad_TrackRes_Completeness
Definition: MuonTrackingEff_module.cc:252
recob::Track::Vertex
Point_t const & Vertex() const
Access to track position at different points.
Definition: Track.h:158
DUNE::MuonTrackingEff::processEff
void processEff(const art::Event &evt, bool &isFiducial)
Definition: MuonTrackingEff_module.cc:1108
sim::ParticleList::begin
iterator begin()
Definition: ParticleList.h:305
Track.h
Provides recob::Track data product.
TrackID
TrackID
Definition: ClusteringValidation_module.cc:58
DUNE::MuonTrackingEff::FuncDistanceAndAngleBetweenTracks
void FuncDistanceAndAngleBetweenTracks(art::Ptr< recob::Track > Track1, art::Ptr< recob::Track > Track2, double &TempDistanceBetweenTracks, double &TempAngleBetweenTracks, double &TempCriteriaTwoTracks)
Definition: MuonTrackingEff_module.cc:1647
geo::TPCID
The data type to uniquely identify a TPC.
Definition: geo_types.h:381
geo_types.h
Definition of data types for geometry description.
hit
Detector simulation of raw signals on wires.
Definition: DumpHits_module.cc:25
DUNE::MuonTrackingEff::h_MuonTrackStitching_FailedCriteria_TrackResSecondMuon_Angle
TH2D * h_MuonTrackStitching_FailedCriteria_TrackResSecondMuon_Angle
Definition: MuonTrackingEff_module.cc:235
cheat::ParticleInventoryService::ParticleList
const sim::ParticleList & ParticleList() const
Definition: ParticleInventoryService.cc:138
simb::MCParticle::Vx
double Vx(const int i=0) const
Definition: MCParticle.h:222
recob::Track::ID
int ID() const
Definition: Track.h:244
DUNE::MuonTrackingEff::h_MuonTrackStitching_FailedCriteriaButLeadingPlusSecondGood_TrackResLeading_TrackResSecond
TH2D * h_MuonTrackStitching_FailedCriteriaButLeadingPlusSecondGood_TrackResLeading_TrackResSecond
Definition: MuonTrackingEff_module.cc:248
art::ProductRetriever::getByLabel
bool getByLabel(std::string const &label, std::string const &instance, Handle< PROD > &result) const
Definition: ProductRetriever.h:444
DUNE::MuonTrackingEff::CountBadLeadingMuonTrackAndLeadingPlusSecondBadTrackTooLong
int CountBadLeadingMuonTrackAndLeadingPlusSecondBadTrackTooLong
Definition: MuonTrackingEff_module.cc:129
detinfo::DetectorClocksData
Contains all timing reference information for the detector.
Definition: DetectorClocksData.h:282
DUNE::MuonTrackingEff::h_MuonTrackStitching_FailedCriteriaAndLeadingPlusSecondBad_Distance_Angle
TH2D * h_MuonTrackStitching_FailedCriteriaAndLeadingPlusSecondBad_Distance_Angle
Definition: MuonTrackingEff_module.cc:254
recob::Track::EndDirection
Vector_t EndDirection() const
Access to track direction at different points.
Definition: Track.h:167
simb::MCParticle::Momentum
const TLorentzVector & Momentum(const int i=0) const
Definition: MCParticle.h:221
DUNE::MuonTrackingEff::h_MuonTrackStitching_CompletenessSecondMuon_Angle
TH2D * h_MuonTrackStitching_CompletenessSecondMuon_Angle
Definition: MuonTrackingEff_module.cc:227
MF_LOG_DEBUG
#define MF_LOG_DEBUG(id)
Definition: MessageLogger.h:294
DUNE::MuonTrackingEff::analyze
void analyze(const art::Event &evt) override
Definition: MuonTrackingEff_module.cc:1100
simb::MCParticle::Vz
double Vz(const int i=0) const
Definition: MCParticle.h:224
recob::Track::End
Point_t const & End() const
Access to track position at different points.
Definition: Track.h:159
DUNE::MuonTrackingEff::h_MuonTrackStitching_FailedCriteriaButLeadingPlusSecondGood_TrackRes_Completeness
TH2D * h_MuonTrackStitching_FailedCriteriaButLeadingPlusSecondGood_TrackRes_Completeness
Definition: MuonTrackingEff_module.cc:246
DUNE::MuonTrackingEff::geom
art::ServiceHandle< geo::Geometry const > geom
Definition: MuonTrackingEff_module.cc:267
DUNE::MuonTrackingEff::truthLength
double truthLength(const simb::MCParticle *MCparticle)
Definition: MuonTrackingEff_module.cc:1709
recob::Hit
2D representation of charge deposited in the TDC/wire plane
Definition: Hit.h:46
DUNE::MuonTrackingEff::h_MuonTrackStitching_FailedCriteria_TrackResLeading_TrackResSecond
TH2D * h_MuonTrackStitching_FailedCriteria_TrackResLeading_TrackResSecond
Definition: MuonTrackingEff_module.cc:232
tca::evt
TCEvent evt
Definition: DataStructs.cxx:8
art::fill_ptr_vector
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
Definition: Ptr.h:306
DUNE::MuonTrackingEff::CountBadLeadingMuonTrackAndOnlyOneMuonTrackTrackTooShort
int CountBadLeadingMuonTrackAndOnlyOneMuonTrackTrackTooShort
Definition: MuonTrackingEff_module.cc:133
geo
Namespace collecting geometry-related classes utilities.
Definition: BeamFlashTrackMatchTaggerAlg.h:20
track
Float_t track
Definition: plot.C:35
Track
Definition: Track.hh:42
DUNE::MuonTrackingEff::CountBadLeadingMuonTrackAndLeadingPlusSecondBadCompleteness
int CountBadLeadingMuonTrackAndLeadingPlusSecondBadCompleteness
Definition: MuonTrackingEff_module.cc:126
simb::MCParticle::Vy
double Vy(const int i=0) const
Definition: MCParticle.h:223
Geometry.h
art framework interface to geometry description
DUNE::MuonTrackingEff::h_MuonTrackStitching_MatchedCriteria_Distance_Angle
TH2D * h_MuonTrackStitching_MatchedCriteria_Distance_Angle
Definition: MuonTrackingEff_module.cc:242
event
Event finding and building.
Definition: EventCheater_module.cc:32
TPCGeo.h
Encapsulate the construction of a single detector plane.
geo::CryostatID
The data type to uniquely identify a cryostat.
Definition: geo_types.h:192
DUNE::MuonTrackingEff::h_MuonTrackStitching_FailedCriteriaButLeadingPlusSecondGood_Distance_Angle
TH2D * h_MuonTrackStitching_FailedCriteriaButLeadingPlusSecondGood_Distance_Angle
Definition: MuonTrackingEff_module.cc:249
vertex
vertex reconstruction
Definition: AggregateVertexAna_module.cc:32