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