LArSoft  v06_85_00
Liquid Argon Software toolkit - http://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 #ifndef MuonTrackingEff_Module
17 #define MuonTrackingEff_Module
18 
19 // LArSoft includes
20 #include "lardata/DetectorInfoServices/DetectorPropertiesService.h"
21 #include "lardata/Utilities/GeometryUtilities.h"
22 #include "larsim/Simulation/LArG4Parameters.h"
23 #include "lardataobj/RecoBase/Track.h"
24 #include "nusimdata/SimulationBase/MCParticle.h"
25 #include "nusimdata/SimulationBase/MCTruth.h"
26 #include "larcoreobj/SimpleTypesAndConstants/geo_types.h"
27 #include "larsim/MCCheater/BackTrackerService.h"
28 #include "larsim/MCCheater/ParticleInventoryService.h"
29 #include "lardata/DetectorInfoServices/DetectorPropertiesService.h"
30 #include "lardata/DetectorInfoServices/DetectorClocksService.h"
31 #include "lardataobj/RawData/ExternalTrigger.h"
32 
33 // Framework includes
34 #include "art/Framework/Core/EDAnalyzer.h"
35 #include "art/Framework/Principal/Event.h"
36 #include "art/Framework/Principal/Handle.h"
37 #include "art/Framework/Services/Registry/ServiceHandle.h"
38 #include "art/Framework/Services/Optional/TFileService.h"
39 #include "art/Framework/Core/ModuleMacros.h"
40 #include "canvas/Persistency/Common/FindManyP.h"
41 #include "canvas/Persistency/Common/PtrVector.h"
42 #include "messagefacility/MessageLogger/MessageLogger.h"
43 #include "fhiclcpp/ParameterSet.h"
44 
45 // ROOT includes
46 #include "TFile.h"
47 #include "TTree.h"
48 #include "TDirectory.h"
49 #include "TH1.h"
50 #include "TH2.h"
51 #include "TEfficiency.h"
52 #include "TGraphAsymmErrors.h"
53 #include "TStyle.h"
54 #include "TVector3.h"
55 
56 #define MAX_TRACKS 1000
57 using namespace std;
58 
59 //========================================================================
60 
61 namespace DUNE{
62 
63 class MuonTrackingEff : public art::EDAnalyzer {
64 public:
65 
66  explicit MuonTrackingEff(fhicl::ParameterSet const& pset);
67  virtual ~MuonTrackingEff();
68 
69  void beginJob();
70  void endJob();
71  void beginRun(const art::Run& run);
72  void analyze(const art::Event& evt);
73 
74  void reconfigure(fhicl::ParameterSet const& pset);
75 
76  void processEff(const art::Event& evt, bool &isFiducial);
77 
78  void truthMatcher( 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);
79 
80  void FuncDistanceAndAngleBetweenTracks(art::Ptr<recob::Track> Track1, art::Ptr<recob::Track> Track2, double &TempDistanceBetweenTracks, double &TempAngleBetweenTracks, double &TempCriteriaTwoTracks);
81 
82  void FuncDistanceAndAngleBetweenTruthAndRecoTrack(const simb::MCParticle *&MCparticle, art::Ptr<recob::Track> Track, double &TempDistanceBetweenTruthAndRecoTrack, double &TempAngleBeetweenTruthAndRecoTrack);
83 
84  double truthLength(const simb::MCParticle *MCparticle);
85 
86  bool insideFV(double vertex[4]);
87 
88  void doEfficiencies();
89 
90 private:
91 
92  // the parameters we'll read from the .fcl
93  std::string fMCTruthModuleLabel;
94  std::string fTrackModuleLabel;
95  int fMuonPDGCode;
96 
97  // int MC_isCC;
98  // int MC_incoming_PDG;
99  // double MC_incoming_P[4];
100  double MCTruthMuonVertex[4];
101  // double MCTruthMuonStartMomentum[4];
102 
103  int MCTruthMuonID;
104  double MCTruthMuonMomentum;
105 
106  // double MCTruthMuonMomentumXZ=0;
107  // double MCTruthMuonMomentumYZ=0;
108  double MCTruthMuonThetaXZ=0;
109  double MCTruthMuonThetaYZ=0;
110 
111  //Counts to calculate efficiency and failed criteria
112  int EventCounter=0;
113 
114  int CountMCTruthMuon=0;
115  int CountRecoMuon=0;
116 
117  int CountGoodLeadingMuonTrack=0;
118  int CountNoRecoTracks=0;
119  int CountNoMuonTracks=0;
120  int CountBadLeadingMuonTrack=0;
121  int CountCompleteness=0;
122  int CountPurity=0;
123  int CountTrackLengthTooShort=0;
124  int CountTrackLengthTooLong=0;
125 
126  int CountBadLeadingMuonTrackAndOnlyOneMuonTrack=0;
127  int CountBadLeadingMuonTrackButLeadingPlusSecondGood=0;
128  int CountBadLeadingMuonTrackAndLeadingPlusSecondBad=0;
129 
130  int CountBadLeadingMuonTrackAndLeadingPlusSecondBadCompleteness=0;
131  int CountBadLeadingMuonTrackAndLeadingPlusSecondBadPurity=0;
132  int CountBadLeadingMuonTrackAndLeadingPlusSecondBadTrackTooShort=0;
133  int CountBadLeadingMuonTrackAndLeadingPlusSecondBadTrackTooLong=0;
134 
135  int CountBadLeadingMuonTrackAndOnlyOneMuonTrackCompleteness=0;
136  int CountBadLeadingMuonTrackAndOnlyOneMuonTrackPurity=0;
137  int CountBadLeadingMuonTrackAndOnlyOneMuonTrackTrackTooShort=0;
138  int CountBadLeadingMuonTrackAndOnlyOneMuonTrackTrackTooLong=0;
139 
140  double Criteria;
141  // int NMuonTracksTooShort=0;
142 
143  int GoodEvents1MuonTrack=0;
144  int GoodEvents2MuonTrack=0;
145  int GoodEvents3MuonTrack=0;
146  int GoodEvents4OrMoreMuonTrack=0;
147 
148  int BadEvents0MuonTrack=0;
149  int BadEvents1MuonTrack=0;
150  int BadEvents2MuonTrack=0;
151  int BadEvents3MuonTrack=0;
152  int BadEvents4OrMoreMuonTrack=0;
153 
154  //TH1Ds
155  //Single Criteria and total reco energy
156  TH1D *h_Purity;
157  TH1D *h_Completeness;
158  TH1D *h_TrackRes;
159  TH1D *h_TotalRecoEnergy;
160  TH1D *h_TruthLength;
161  TH1D *h_VertexRes;
162  TH1D *h_DirectionRes;
163 
164  //Efficiency ThetaXZ
165  TH1D *h_Efficiency_ThetaXZ;
166  TH1D *h_ThetaXZ_den;
167  TH1D *h_ThetaXZ_num;
168 
169  //Efficiency ThetaYZ
170  TH1D *h_Efficiency_ThetaYZ;
171  TH1D *h_ThetaYZ_den;
172  TH1D *h_ThetaYZ_num;
173 
174  //Efficiency SinThetaYZ
175  TH1D *h_Efficiency_SinThetaYZ;
176  TH1D *h_SinThetaYZ_den;
177  TH1D *h_SinThetaYZ_num;
178 
179 
180  //TH2Ds
181  //Efficiency ThetaXZ vs. ThetaYZ
182  TH2D *h_Efficiency_ThetaXZ_ThetaYZ;
183  TH2D *h_ThetaXZ_ThetaYZ_den;
184  TH2D *h_ThetaXZ_ThetaYZ_num;
185  TH2D *h_FailedReconstruction_ThetaXZ_ThetaYZ;
186 
187  //Efficiency ThetaXZ vs. SinThetaYZ
188  TH2D *h_Efficiency_ThetaXZ_SinThetaYZ;
189  TH2D *h_ThetaXZ_SinThetaYZ_den;
190  TH2D *h_ThetaXZ_SinThetaYZ_num;
191  TH2D *h_FailedReconstruction_ThetaXZ_SinThetaYZ;
192 
193  //Efficiency ThetaXZ vs. ThetaYZ after summing up leading plus second
194  TH2D *h_Efficiency_ThetaXZ_ThetaYZ_LeadingPlusSecond;
195  TH2D *h_ThetaXZ_ThetaYZ_LeadingPlusSecondOk;
196  TH2D *h_ThetaXZ_ThetaYZ_LeadingPlusSecondOk_num;
197 
198  //Difference in efficiency before and after summing up leading plus second: ThetaXZ vs. ThetaYZ
199  TH2D *h_Efficiency_ThetaXZ_ThetaYZ_DifferenceLeadingAndLeadingPlusSecond;
200 
201  //Efficiency ThetaXZ vs. SinThetaYZ after summing up leading plus second
202  TH2D *h_Efficiency_ThetaXZ_SinThetaYZ_LeadingPlusSecond;
203  TH2D *h_ThetaXZ_SinThetaYZ_LeadingPlusSecondOk;
204  TH2D *h_ThetaXZ_SinThetaYZ_LeadingPlusSecondOk_num;
205 
206  //Failed Criteria
207  TH2D *h_NoRecoTrackAtAll_ThetaXZ_SinThetaYZ;
208  TH2D *h_NoMuonTrack_ThetaXZ_SinThetaYZ;
209  TH2D *h_TrackTooShort_ThetaXZ_SinThetaYZ;
210  TH2D *h_TrackTooLong_ThetaXZ_SinThetaYZ;
211  TH2D *h_Completeness_ThetaXZ_SinThetaYZ;
212  TH2D *h_Purity_ThetaXZ_SinThetaYZ;
213 
214  //Criteria vs. NRecoTrack
215  TH2D *h_Criteria_NRecoTrack;
216  TH2D *h_Criteria_NRecoTrack_den;
217  TH2D *h_Criteria_NRecoTrack_num;
218 
219  //Criteria vs. NMuonTrack
220  TH2D *h_Criteria_NMuonTrack;
221  TH2D *h_Criteria_NMuonTrack_den;
222  TH2D *h_Criteria_NMuonTrack_num;
223 
224  //NoMuonTrack: Max length of no muon track vs. PDG code of that track (MC truth)
225  TH2D *h_NoMuonTrack_MaxTrackLength_PDGCode;
226 
227  //Stitching variables: all events
228  TH2D *h_MuonTrackStitching_TrackRes_Completeness;
229  TH2D *h_MuonTrackStitching_TrackResLeading_TrackResSecond;
230  TH2D *h_MuonTrackStitching_Distance_Angle;
231  TH2D *h_MuonTrackStitching_TrackResSecondMuon_Angle;
232  TH2D *h_MuonTrackStitching_CompletenessSecondMuon_Angle;
233  TH2D *h_MuonTrackStitching_CriteriaTwoTracks_Angle;
234 
235  //Stitching variables: bad events
236  TH2D *h_MuonTrackStitching_FailedCriteria_TrackRes_Completeness;
237  TH2D *h_MuonTrackStitching_FailedCriteria_TrackResLeading_TrackResSecond;
238  TH2D *h_MuonTrackStitching_FailedCriteria_CompletenessLeading_CompletenessSecond;
239  TH2D *h_MuonTrackStitching_FailedCriteria_Distance_Angle;
240  TH2D *h_MuonTrackStitching_FailedCriteria_TrackResSecondMuon_Angle;
241  TH2D *h_MuonTrackStitching_FailedCriteria_CompletenessSecondMuon_Angle;
242  TH2D *h_MuonTrackStitching_FailedCriteria_CriteriaTwoTracks_Angle;
243 
244  //Stitching variables: good events
245  TH2D *h_MuonTrackStitching_MatchedCriteria_TrackRes_Completeness;
246  TH2D *h_MuonTrackStitching_MatchedCriteria_TrackResLeading_TrackResSecond;
247  TH2D *h_MuonTrackStitching_MatchedCriteria_Distance_Angle;
248  TH2D *h_MuonTrackStitching_MatchedCriteria_CriteriaTwoTracks_Angle;
249 
250  //Stitching variables: bad events but leading plus second ok
251  TH2D *h_MuonTrackStitching_FailedCriteriaButLeadingPlusSecondGood_TrackRes_Completeness;
252  TH2D *h_MuonTrackStitching_FailedCriteriaButLeadingPlusSecondGood_TrackResLeading_TrackResSecond;
253  TH2D *h_MuonTrackStitching_FailedCriteriaButLeadingPlusSecondGood_Distance_Angle;
254 
255  //Stitching variables: bad events and leading plus second not ok
256  TH2D *h_MuonTrackStitching_FailedCriteriaAndLeadingPlusSecondBad_TrackRes_Completeness;
257  TH2D *h_MuonTrackStitching_FailedCriteriaAndLeadingPlusSecondBad_TrackResLeading_TrackResSecond;
258  TH2D *h_MuonTrackStitching_FailedCriteriaAndLeadingPlusSecondBad_Distance_Angle;
259 
260 
261  float fFidVolCutX;
262  float fFidVolCutY;
263  float fFidVolCutZ;
264 
265  float fFidVolXmin;
266  float fFidVolXmax;
267  float fFidVolYmin;
268  float fFidVolYmax;
269  float fFidVolZmin;
270  float fFidVolZmax;
271 
272  detinfo::DetectorProperties const *detprop = lar::providerFrom<detinfo::DetectorPropertiesService>();
273  detinfo::DetectorClocks const *ts = lar::providerFrom<detinfo::DetectorClocksService>();
274  double XDriftVelocity = detprop->DriftVelocity()*1e-3; //cm/ns
275  double WindowSize = detprop->NumberTimeSamples() * ts->TPCClock().TickPeriod() * 1e3;
276  art::ServiceHandle<geo::Geometry> geom;
277 
278  //My histograms
279  int NThetaXZBins=36;
280  int ThetaXZBinMin=0;
281  int ThetaXZBinMax=360;
282 
283  int NThetaYZBins=18;
284  int ThetaYZBinMin=-90;
285  int ThetaYZBinMax=90;
286 
287  int NSinThetaYZBins=18;
288  int SinThetaYZBinMin=-1;
289  int SinThetaYZBinMax=1;
290 
291  int NCriteriaBins=13;
292  double CriteriaBinMin=-0.25;
293  double CriteriaBinMax=6.25;
294 
295  int NRecoTracksBins=19;
296  double RecoTracksBinMin=-0.25;
297  double RecoTracksBinMax=9.25;
298 
299  }; // class MuonTrackingEff
300 
301 
302 //========================================================================
303 MuonTrackingEff::MuonTrackingEff(fhicl::ParameterSet const& parameterSet)
304  : EDAnalyzer(parameterSet)
305 {
306  reconfigure(parameterSet);
307 }
308 //========================================================================
309 MuonTrackingEff::~MuonTrackingEff(){
310  //destructor
311 }
312 //========================================================================
313 void MuonTrackingEff::reconfigure(fhicl::ParameterSet const& p){
314 
315  fMCTruthModuleLabel = p.get<std::string>("MCTruthModuleLabel");
316  fTrackModuleLabel = p.get<std::string>("TrackModuleLabel");
317  fMuonPDGCode = p.get<int>("MuonPDGCode");
318  fFidVolCutX = p.get<float>("FidVolCutX");
319  fFidVolCutY = p.get<float>("FidVolCutY");
320  fFidVolCutZ = p.get<float>("FidVolCutZ");
321 }
322 //========================================================================
323 void MuonTrackingEff::beginJob(){
324  std::cout<<"job begin..."<<std::endl;
325  // Get geometry.
326  auto const* geo = lar::providerFrom<geo::Geometry>();
327  // Define histogram boundaries (cm).
328  // For now only draw cryostat=0.
329  double minx = 1e9;
330  double maxx = -1e9;
331  double miny = 1e9;
332  double maxy = -1e9;
333  double minz = 1e9;
334  double maxz = -1e9;
335  for (size_t i = 0; i<geo->NTPC(); ++i){
336  double local[3] = {0.,0.,0.};
337  double world[3] = {0.,0.,0.};
338  const geo::TPCGeo &tpc = geo->TPC(i);
339  tpc.LocalToWorld(local,world);
340  if (minx>world[0]-geo->DetHalfWidth(i))
341  minx = world[0]-geo->DetHalfWidth(i);
342  if (maxx<world[0]+geo->DetHalfWidth(i))
343  maxx = world[0]+geo->DetHalfWidth(i);
344  if (miny>world[1]-geo->DetHalfHeight(i))
345  miny = world[1]-geo->DetHalfHeight(i);
346  if (maxy<world[1]+geo->DetHalfHeight(i))
347  maxy = world[1]+geo->DetHalfHeight(i);
348  if (minz>world[2]-geo->DetLength(i)/2.)
349  minz = world[2]-geo->DetLength(i)/2.;
350  if (maxz<world[2]+geo->DetLength(i)/2.)
351  maxz = world[2]+geo->DetLength(i)/2.;
352  }
353 
354  fFidVolXmin = minx + fFidVolCutX;
355  fFidVolXmax = maxx - fFidVolCutX;
356  fFidVolYmin = miny + fFidVolCutY;
357  fFidVolYmax = maxy - fFidVolCutY;
358  fFidVolZmin = minz + fFidVolCutZ;
359  fFidVolZmax = maxz - fFidVolCutZ;
360 
361  std::cout<<"Fiducial volume:"<<"\n"
362  <<fFidVolXmin<<"\t< x <\t"<<fFidVolXmax<<"\n"
363  <<fFidVolYmin<<"\t< y <\t"<<fFidVolYmax<<"\n"
364  <<fFidVolZmin<<"\t< z <\t"<<fFidVolZmax<<"\n";
365 
366  art::ServiceHandle<art::TFileService> tfs;
367 
368  //TH1D's
369  gStyle->SetTitleOffset (1.3,"Y");
370 
371  //Single Criteria and total reco energy
372  h_Purity = tfs->make<TH1D>("h_Purity","All events: Purity vs. # events; Purity; # events",60,0,1.2);
373 
374  h_Completeness = tfs->make<TH1D>("h_Completeness","All events: Completeness vs # events; Completeness; # events",60,0,1.2);
375  h_Completeness->SetLineColor(kBlue);
376 
377  h_TrackRes = tfs->make<TH1D>("h_TrackRes", "All events: L_{reco}/L_{truth} vs. # events; L_{reco}/L_{truth}; # events;",75,0,1.5);
378  h_TrackRes->SetLineColor(kRed);
379 
380  h_TotalRecoEnergy = tfs->make<TH1D>("h_TotalRecoEnergy", "All events: Total reco energy (sum of all hits in all tracks) vs. # events; E_{reco., tot.} [MeV]; # events",100,0,1000);
381 
382  h_TruthLength = tfs->make<TH1D>("h_TruthLength", "All events: truth muon length vs. # events; truth muon length [cm]; # events",100,0,300);
383 
384  h_VertexRes = tfs->make<TH1D>("h_VertexRes", "All events: Vertex residuals vs. # events; #Delta vertex_{truth-teco} [cm]; # events", 300,0,300);
385 
386  h_DirectionRes = tfs->make<TH1D>("h_DirectionRes", "All events: Angular residuals vs. # events; #Delta#theta_{truth-reco} [#circ]; # events", 180,0,180);
387 
388  //Efficiency ThetaXZ
389  h_Efficiency_ThetaXZ = tfs->make<TH1D>("h_Efficiency_ThetaXZ","Muon reco efficiency vs. #theta_{XZ}; #theta_{XZ} [#circ]; Efficiency",NThetaXZBins,ThetaXZBinMin,ThetaXZBinMax);
390  h_ThetaXZ_den = tfs->make<TH1D>("h_ThetaXZ_den","# generated muons vs. #theta_{XZ}; #theta_{XZ} [#circ]; # generated muons",NThetaXZBins,ThetaXZBinMin,ThetaXZBinMax);
391  h_ThetaXZ_num = tfs->make<TH1D>("h_ThetaXZ_num","# reconstructed muons vs. #theta_{XZ}; #theta_{XZ} [#circ]; # reconstructed muons",NThetaXZBins,ThetaXZBinMin,ThetaXZBinMax);
392 
393  //Efficiency ThetaYZ
394  h_Efficiency_ThetaYZ = tfs->make<TH1D>("h_Efficiency_ThetaYZ","Muon reco efficiency vs. #theta_{YZ}; #theta_{YZ} [#circ]; Muon reco efficiency",NThetaYZBins,ThetaYZBinMin,ThetaYZBinMax);;
395  h_ThetaYZ_den = tfs->make<TH1D>("h_ThetaYZ_den","# generated muons vs. #theta_{YZ}; #theta_{YZ} [#circ]; # generated muons",NThetaYZBins,ThetaYZBinMin,ThetaYZBinMax);
396  h_ThetaYZ_num = tfs->make<TH1D>("h_ThetaYZ_num","# reconstructed muons vs. #theta_{YZ}; #theta_{YZ} [#circ]; # reconstructed muons",NThetaYZBins,ThetaYZBinMin,ThetaYZBinMax);
397 
398  //Efficiency SinThetaYZ
399  h_Efficiency_SinThetaYZ = tfs->make<TH1D>("h_Efficiency_SinThetaYZ","Muon reco efficiency vs. sin(#theta_{YZ}); sin(#theta_{YZ}); Muon reco efficiency",NSinThetaYZBins,SinThetaYZBinMin,SinThetaYZBinMax);;
400  h_SinThetaYZ_den = tfs->make<TH1D>("h_SinThetaYZ_den","# generated muons vs. sin(#theta_{YZ}); sin(#theta_{YZ}); # generated muons",NSinThetaYZBins,SinThetaYZBinMin,SinThetaYZBinMax);
401  h_SinThetaYZ_num = tfs->make<TH1D>("h_SinThetaYZ_num","# reconstructed muons vs. sin(#theta_{YZ}); sin(#theta_{YZ}); # reconstructed muons",NSinThetaYZBins,SinThetaYZBinMin,SinThetaYZBinMax);
402 
403  h_Purity->Sumw2();
404  h_Completeness->Sumw2();
405  h_TrackRes->Sumw2();
406  h_TotalRecoEnergy->Sumw2();
407  h_TruthLength->Sumw2();
408  h_VertexRes->Sumw2();
409  h_DirectionRes->Sumw2();
410 
411  h_Efficiency_SinThetaYZ->Sumw2();
412  h_SinThetaYZ_den->Sumw2();
413  h_SinThetaYZ_num->Sumw2();
414 
415  h_Efficiency_ThetaXZ->Sumw2();
416  h_ThetaXZ_den->Sumw2();
417  h_ThetaXZ_num->Sumw2();
418 
419  h_Efficiency_ThetaYZ->Sumw2();
420  h_ThetaYZ_den->Sumw2();
421  h_ThetaYZ_num->Sumw2();
422 
423 
424  //TH2D's
425  //Efficiency and Failed Reconstruction ThetaXZ vs. ThetaYZ
426  h_Efficiency_ThetaXZ_ThetaYZ = tfs->make<TH2D>("h_Efficiency_ThetaXZ_ThetaYZ","Muon reco efficiency: #theta_{XZ} vs. #theta_{YZ}; #theta_{XZ} [#circ]; #theta_{YZ} [#circ]",NThetaXZBins,ThetaXZBinMin,ThetaXZBinMax,NThetaYZBins,ThetaYZBinMin,ThetaYZBinMax);
427  h_Efficiency_ThetaXZ_ThetaYZ->SetOption("colz");
428 
429  h_ThetaXZ_ThetaYZ_den = tfs->make<TH2D>("h_ThetaXZ_ThetaYZ_den","# generated muons: #theta_{XZ} vs. #theta_{YZ}; #theta_{XZ} [#circ]; #theta_{YZ} [#circ]",NThetaXZBins,ThetaXZBinMin,ThetaXZBinMax,NThetaYZBins,ThetaYZBinMin,ThetaYZBinMax);
430  h_ThetaXZ_ThetaYZ_den->SetOption("colz");
431 
432  h_ThetaXZ_ThetaYZ_num = tfs->make<TH2D>("h_ThetaXZ_ThetaYZ_num","# reconstructed muons: #theta_{XZ} vs. #theta_{YZ}; #theta_{XZ} [#circ]; #theta_{YZ} [#circ]",NThetaXZBins,ThetaXZBinMin,ThetaXZBinMax,NThetaYZBins,ThetaYZBinMin,ThetaYZBinMax);
433  h_ThetaXZ_ThetaYZ_num->SetOption("colz");
434 
435  h_FailedReconstruction_ThetaXZ_ThetaYZ = tfs->make<TH2D>("h_FailedReconstruction_ThetaXZ_ThetaYZ","# failed reconstructions: #theta_{XZ} vs. #theta_{YZ}; #theta_{XZ} [#circ]; #theta_{YZ} [#circ]",NThetaXZBins,ThetaXZBinMin,ThetaXZBinMax,NThetaYZBins,ThetaYZBinMin,ThetaYZBinMax);
436  h_FailedReconstruction_ThetaXZ_ThetaYZ->SetOption("colz");
437 
438  //Efficiency and Failed Reconstruction ThetaXZ vs. SinThetaYZ
439  h_Efficiency_ThetaXZ_SinThetaYZ = tfs->make<TH2D>("h_Efficiency_ThetaXZ_SinThetaYZ","Muon reco efficiency: #theta_{XZ} vs. sin(#theta_{YZ}); #theta_{XZ} [#circ]; sin(#theta_{YZ})",NThetaXZBins,ThetaXZBinMin,ThetaXZBinMax,NSinThetaYZBins,SinThetaYZBinMin,SinThetaYZBinMax);
440  h_Efficiency_ThetaXZ_SinThetaYZ->SetOption("colz");
441 
442  h_ThetaXZ_SinThetaYZ_den = tfs->make<TH2D>("h_ThetaXZ_SinThetaYZ_den","# generated muons: #theta_{XZ} vs. sin(#theta_{YZ}); #theta_{XZ} [#circ]; sin(#theta_{YZ})",NThetaXZBins,ThetaXZBinMin,ThetaXZBinMax,NSinThetaYZBins,SinThetaYZBinMin,SinThetaYZBinMax);
443  h_ThetaXZ_SinThetaYZ_den->SetOption("colz");
444 
445  h_ThetaXZ_SinThetaYZ_num = tfs->make<TH2D>("h_ThetaXZ_SinThetaYZ_num","# reconstructed muons: #theta_{XZ} vs. sin(#theta_{YZ}); #theta_{XZ} [#circ]; sin(#theta_{YZ}); #theta_{XZ} [#circ]; sin(#theta_{YZ})",NThetaXZBins,ThetaXZBinMin,ThetaXZBinMax,NSinThetaYZBins,SinThetaYZBinMin,SinThetaYZBinMax);
446  h_ThetaXZ_SinThetaYZ_num->SetOption("colz");
447 
448  h_FailedReconstruction_ThetaXZ_SinThetaYZ = tfs->make<TH2D>("h_FailedReconstruction_ThetaXZ_SinThetaYZ","# failed reconstructions: #theta_{XZ} vs. sin(#theta_{YZ}); #theta_{XZ} [#circ]; sin(#theta_{YZ})",NThetaXZBins,ThetaXZBinMin,ThetaXZBinMax,NSinThetaYZBins,SinThetaYZBinMin,SinThetaYZBinMax);
449  h_FailedReconstruction_ThetaXZ_SinThetaYZ->SetOption("colz");
450 
451  //Efficiency ThetaXZ vs. ThetaYZ after summing up leading plus second
452  h_Efficiency_ThetaXZ_ThetaYZ_LeadingPlusSecond = tfs->make<TH2D>("h_Efficiency_ThetaXZ_ThetaYZ_LeadingPlusSecond","Muon reco efficiency after stitching: #theta_{XZ} vs. #theta_{YZ}; #theta_{XZ} [#circ]; #theta_{YZ} [#circ]",NThetaXZBins,ThetaXZBinMin,ThetaXZBinMax,NThetaYZBins,ThetaYZBinMin,ThetaYZBinMax);
453  h_Efficiency_ThetaXZ_ThetaYZ_LeadingPlusSecond->SetOption("colz");
454 
455  h_ThetaXZ_ThetaYZ_LeadingPlusSecondOk = tfs->make<TH2D>("h_ThetaXZ_ThetaYZ_LeadingPlusSecondOk","# reconstructed muons after stitching (failed before stitching): #theta_{XZ} vs #theta_{YZ}; #theta_{XZ} [#circ]; #theta_{YZ} [#circ]",NThetaXZBins,ThetaXZBinMin,ThetaXZBinMax,NThetaYZBins,ThetaYZBinMin,ThetaYZBinMax);
456  h_ThetaXZ_ThetaYZ_LeadingPlusSecondOk->SetOption("colz");
457 
458  h_ThetaXZ_ThetaYZ_LeadingPlusSecondOk_num = tfs->make<TH2D>("h_ThetaXZ_ThetaYZ_LeadingPlusSecondOk_num","# reconstructed muons after stitching: #theta_{XZ} vs. #theta_{YZ}; #theta_{XZ} [#circ]; #theta_{YZ} [#circ]",NThetaXZBins,ThetaXZBinMin,ThetaXZBinMax,NThetaYZBins,ThetaYZBinMin,ThetaYZBinMax);
459  h_ThetaXZ_ThetaYZ_LeadingPlusSecondOk_num->SetOption("colz");
460 
461  //Efficiency ThetaXZ vs. SinThetaYZ after summing up leading plus second
462  h_Efficiency_ThetaXZ_SinThetaYZ_LeadingPlusSecond = tfs->make<TH2D>("h_Efficiency_ThetaXZ_SinThetaYZ_LeadingPlusSecond","Muon reco efficiency after stitching: #theta_{XZ} vs. sin(#theta_{YZ}); #theta_{XZ} [#circ]; sin(#theta_{YZ})",NThetaXZBins,ThetaXZBinMin,ThetaXZBinMax,NSinThetaYZBins,SinThetaYZBinMin,SinThetaYZBinMax);
463  h_Efficiency_ThetaXZ_SinThetaYZ_LeadingPlusSecond->SetOption("colz");
464 
465  h_ThetaXZ_SinThetaYZ_LeadingPlusSecondOk = tfs->make<TH2D>("h_ThetaXZ_SinThetaYZ_LeadingPlusSecondOk","# reconstructed muons after stitching (failed before stitching): #theta_{XZ} vs. sin(#theta_{YZ}); #theta_{XZ} [#circ]; sin(#theta_{YZ})",NThetaXZBins,ThetaXZBinMin,ThetaXZBinMax,NSinThetaYZBins,SinThetaYZBinMin,SinThetaYZBinMax);
466  h_ThetaXZ_SinThetaYZ_LeadingPlusSecondOk->SetOption("colz");
467 
468  h_ThetaXZ_SinThetaYZ_LeadingPlusSecondOk_num = tfs->make<TH2D>("h_ThetaXZ_SinThetaYZ_LeadingPlusSecondOk_num","# reconstructed muons after stitching: #theta_{XZ} vs. sin(#theta_{YZ}); #theta_{XZ} [#circ]; sin(#theta_{YZ})",NThetaXZBins,ThetaXZBinMin,ThetaXZBinMax,NSinThetaYZBins,SinThetaYZBinMin,SinThetaYZBinMax);
469  h_ThetaXZ_SinThetaYZ_LeadingPlusSecondOk_num->SetOption("colz");
470 
471  //Difference in efficiency before and after summing up leading plus second: ThetaXZ vs. ThetaYZ
472  h_Efficiency_ThetaXZ_ThetaYZ_DifferenceLeadingAndLeadingPlusSecond = tfs->make<TH2D>("h_Efficiency_ThetaXZ_ThetaYZ_DifferenceLeadingAndLeadingPlusSecond","Muon reco efficiency: difference before and after stitching: #theta_{XZ} vs. #theta_{YZ}; #theta_{XZ} [#circ]; #theta_{YZ} [#circ]",NThetaXZBins,ThetaXZBinMin,ThetaXZBinMax,NThetaYZBins,ThetaYZBinMin,ThetaYZBinMax);
473  h_Efficiency_ThetaXZ_ThetaYZ_DifferenceLeadingAndLeadingPlusSecond->SetOption("colz");
474 
475  //Failed Criteria
476  h_NoRecoTrackAtAll_ThetaXZ_SinThetaYZ = tfs->make<TH2D>("h_NoRecoTrackAtAll_ThetaXZ_SinThetaYZ","# events with no reco track at all: #theta_{XZ} vs. sin(#theta_{YZ}); #theta_{XZ} [#circ]; sin(#theta_{YZ})",NThetaXZBins,ThetaXZBinMin,ThetaXZBinMax,NSinThetaYZBins,SinThetaYZBinMin,SinThetaYZBinMax);
477  h_NoRecoTrackAtAll_ThetaXZ_SinThetaYZ->SetOption("colz");
478 
479  h_NoMuonTrack_ThetaXZ_SinThetaYZ = tfs->make<TH2D>("h_NoMuonTrack_ThetaXZ_SinThetaYZ","# events with no muon track: #theta_{XZ} vs. sin(#theta_{YZ}); #theta_{XZ} [#circ]; sin(#theta_{YZ})",NThetaXZBins,ThetaXZBinMin,ThetaXZBinMax,NSinThetaYZBins,SinThetaYZBinMin,SinThetaYZBinMax);
480  h_NoMuonTrack_ThetaXZ_SinThetaYZ->SetOption("colz");
481 
482  h_TrackTooShort_ThetaXZ_SinThetaYZ = tfs->make<TH2D>("h_TrackTooShort_ThetaXZ_SinThetaYZ","# events with L_{reco}/L_{truth} < 75%: #theta_{XZ} vs. sin(#theta_{YZ}); #theta_{XZ} [#circ]; sin(#theta_{YZ})",NThetaXZBins,ThetaXZBinMin,ThetaXZBinMax,NSinThetaYZBins,SinThetaYZBinMin,SinThetaYZBinMax);
483  h_TrackTooShort_ThetaXZ_SinThetaYZ->SetOption("colz");
484 
485  h_TrackTooLong_ThetaXZ_SinThetaYZ = tfs->make<TH2D>("h_TrackTooLong_ThetaXZ_SinThetaYZ","#events with L_{reco}/L_{truth} > 125%: #theta_{XZ} vs. sin(#theta_{YZ}); #theta_{XZ} [#circ]; sin(#theta_{YZ})",NThetaXZBins,ThetaXZBinMin,ThetaXZBinMax,NSinThetaYZBins,SinThetaYZBinMin,SinThetaYZBinMax);
486  h_TrackTooLong_ThetaXZ_SinThetaYZ->SetOption("colz");
487 
488  h_Completeness_ThetaXZ_SinThetaYZ = tfs->make<TH2D>("h_Completeness_ThetaXZ_SinThetaYZ","# events with Completeness < 50%: #theta_{XZ} vs. sin(#theta_{YZ}); #theta_{XZ} [#circ]; sin(#theta_{YZ})",NThetaXZBins,ThetaXZBinMin,ThetaXZBinMax,NSinThetaYZBins,SinThetaYZBinMin,SinThetaYZBinMax);
489  h_Completeness_ThetaXZ_SinThetaYZ->SetOption("colz");
490 
491  h_Purity_ThetaXZ_SinThetaYZ = tfs->make<TH2D>("h_Purity_ThetaXZ_SinThetaYZ","# events with Purity < 50%: #theta_{XZ} vs. sin(#theta_{YZ}); #theta_{XZ} [#circ]; sin(#theta_{YZ})",NThetaXZBins,ThetaXZBinMin,ThetaXZBinMax,NSinThetaYZBins,SinThetaYZBinMin,SinThetaYZBinMax);
492  h_Purity_ThetaXZ_SinThetaYZ->SetOption("colz");
493 
494  //Criteria vs. NRecoTrack
495  h_Criteria_NRecoTrack = tfs->make<TH2D>("h_Criteria_NRecoTrack","Ratio: criteria vs. # reco tracks; Criteria; # reco tracks",NCriteriaBins,CriteriaBinMin,CriteriaBinMax,NRecoTracksBins,RecoTracksBinMin,RecoTracksBinMax);
496  h_Criteria_NRecoTrack->SetOption("colz");
497 
498  h_Criteria_NRecoTrack_num = tfs->make<TH2D>("h_Criteria_NRecoTrack_num","# events: criteria vs. # reco tracks; Criteria; # reco tracks",NCriteriaBins,CriteriaBinMin,CriteriaBinMax,NRecoTracksBins,RecoTracksBinMin,RecoTracksBinMax);
499  h_Criteria_NRecoTrack_num->SetOption("colz");
500 
501  h_Criteria_NRecoTrack_den = tfs->make<TH2D>("h_Criteria_NRecoTrack_den","Divider histogram; Criteria; # reco tracks",NCriteriaBins,CriteriaBinMin,CriteriaBinMax,NRecoTracksBins,RecoTracksBinMin,RecoTracksBinMax);
502  h_Criteria_NRecoTrack_den->SetOption("colz");
503 
504  //Criteria vs. NMuonTrack
505  h_Criteria_NMuonTrack = tfs->make<TH2D>("h_Criteria_NMuonTrack","Ratio: criteria vs. # muon tracks; Criteria; # muon tracks",NCriteriaBins,CriteriaBinMin,CriteriaBinMax,NRecoTracksBins,RecoTracksBinMin,RecoTracksBinMax);
506  h_Criteria_NMuonTrack->SetOption("colz");
507 
508  h_Criteria_NMuonTrack_num = tfs->make<TH2D>("h_Criteria_NMuonTrack_num","# events: criteria vs. # muon tracks; Criteria; # muon tracks",NCriteriaBins,CriteriaBinMin,CriteriaBinMax,NRecoTracksBins,RecoTracksBinMin,RecoTracksBinMax);
509  h_Criteria_NMuonTrack_num->SetOption("colz");
510 
511  h_Criteria_NMuonTrack_den = tfs->make<TH2D>("h_Criteria_NMuonTrack_den","Divider histogram; Criteria; # muon tracks",NCriteriaBins,CriteriaBinMin,CriteriaBinMax,NRecoTracksBins,RecoTracksBinMin,RecoTracksBinMax);
512  h_Criteria_NMuonTrack_den->SetOption("colz");
513 
514  //NoMuonTrack: Max length of no muon track vs. PDG code of that track (MC truth)
515  h_NoMuonTrack_MaxTrackLength_PDGCode = tfs->make<TH2D>("h_NoMuonTrack_MaxTrackLength_PDGCode", "Events with no muon track: L_{reco, max} vs. PDG Code; L_{reco} [cm]; PDG Code",100,0.,200.,200,-100.,100.);
516  h_NoMuonTrack_MaxTrackLength_PDGCode->SetOption("colz");
517 
518  //Stitching variables: all events
519  h_MuonTrackStitching_TrackRes_Completeness = tfs->make<TH2D>("h_MuonTrackStitching_TrackRes_Completeness", "All events: L_{reco}/L_{truth} (leading) vs. Completeness (leading); L_{reco}/L_{truth} (leading); Completeness (leading)",150,0.,1.5,150,0.,1.5);
520  h_MuonTrackStitching_TrackRes_Completeness->SetOption("colz");
521 
522  h_MuonTrackStitching_TrackResLeading_TrackResSecond = tfs->make<TH2D>("h_MuonTrackStitching_TrackResLeading_TrackResSecond", "All events: L_{reco}/L_{truth} (leading) vs. L_{reco}/L_{truth} (second); L_{reco}/L_{truth} (leading); L_{reco}/L_{truth} (second)",150,0.,1.5,150,0.,1.5);
523  h_MuonTrackStitching_TrackResLeading_TrackResSecond->SetOption("colz");
524 
525  h_MuonTrackStitching_Distance_Angle = tfs->make<TH2D>("h_MuonTrackStitching_Distance_Angle", "All events: distance vs. angle b/w leading and second muon track; Distance [cm]; Angle [#circ]",100,0.,100.,100,0.,180.);
526  h_MuonTrackStitching_Distance_Angle->SetOption("colz");
527 
528  h_MuonTrackStitching_TrackResSecondMuon_Angle = tfs->make<TH2D>("h_MuonTrackStitching_TrackResSecondMuon_Angle", "All events: L_{reco}/L_{truth} (second) vs. angle; L_{reco}/L_{truth} (second); Angle [#circ]",150,0.,1.5,180,0,180.);
529  h_MuonTrackStitching_TrackResSecondMuon_Angle->SetOption("colz");
530 
531  h_MuonTrackStitching_CompletenessSecondMuon_Angle = tfs->make<TH2D>("h_MuonTrackStitching_CompletenessSecondMuon_Angle", "All events: Completeness (second) vs. angle; Completeness (second); Angle [#circ]",120,0.,1.2,180,0,180.);
532  h_MuonTrackStitching_CompletenessSecondMuon_Angle->SetOption("colz");
533 
534  h_MuonTrackStitching_CriteriaTwoTracks_Angle = tfs->make<TH2D>("h_MuonTrackStitching_CriteriaTwoTracks_Angle", "All events: CriteriaTwoTracks vs. angle; Criteria; Angle [#circ]",7,0.75,4.25,180,0,180.);
535  h_MuonTrackStitching_CriteriaTwoTracks_Angle->SetOption("colz");
536 
537  //Stitching variables: bad events
538  h_MuonTrackStitching_FailedCriteria_TrackRes_Completeness = tfs->make<TH2D>("h_MuonTrackStitching_FailedCriteria_TrackRes_Completeness", "Bad events: L_{reco}/L_{truth} (leading) vs. Completeness (leading); L_{reco}/L_{truth} (leading); Completeness (leading)",150,0.,1.5,150,0.,1.5);
539  h_MuonTrackStitching_FailedCriteria_TrackRes_Completeness->SetOption("colz");
540 
541  h_MuonTrackStitching_FailedCriteria_Distance_Angle = tfs->make<TH2D>("h_MuonTrackStitching_FailedCriteria_Distance_Angle", "Bad events: distance vs. angle b/w leading and second muon track; Distance [cm]; Angle [#circ]",100,0.,100.,100,0.,180.);
542  h_MuonTrackStitching_FailedCriteria_Distance_Angle->SetOption("colz");
543 
544  h_MuonTrackStitching_FailedCriteria_TrackResLeading_TrackResSecond = tfs->make<TH2D>("h_MuonTrackStitching_FailedCriteria_TrackResLeading_TrackResSecond", "Bad events: L_{reco}/L_{truth} (leading) vs. L_{reco}/L_{truth} (second); L_{reco}/L_{truth} (leading); L_{reco}/L_{truth} (second)",150,0.,1.5,150,0.,1.5);
545  h_MuonTrackStitching_FailedCriteria_TrackResLeading_TrackResSecond->SetOption("colz");
546 
547  h_MuonTrackStitching_FailedCriteria_CompletenessLeading_CompletenessSecond = tfs->make<TH2D>("*h_MuonTrackStitching_FailedCriteria_CompletenessLeading_CompletenessSecond", "Bad events: Completeness (leading) vs. Completeness (second); Completeness (leading); Completeness (second)",150,0.,1.5,150,0.,1.5);
548  h_MuonTrackStitching_FailedCriteria_CompletenessLeading_CompletenessSecond->SetOption("colz");
549 
550  h_MuonTrackStitching_FailedCriteria_TrackResSecondMuon_Angle = tfs->make<TH2D>("h_MuonTrackStitching_FailedCriteria_TrackResSecondMuon_Angle", "Bad events: L_{reco}/L_{truth} (second) vs. angle; L_{reco}/L_{truth} (second); Angle [#circ]",150,0.,1.5,180,0,180.);
551  h_MuonTrackStitching_FailedCriteria_TrackResSecondMuon_Angle->SetOption("colz");
552 
553  h_MuonTrackStitching_FailedCriteria_CompletenessSecondMuon_Angle = tfs->make<TH2D>("h_MuonTrackStitching_FailedCriteria_CompletenessSecondMuon_Angle", "Bad events: Completeness (second) vs. angle; Completeness (second); Angle [#circ]",120,0.,1.2,180,0,180.);
554  h_MuonTrackStitching_FailedCriteria_CompletenessSecondMuon_Angle->SetOption("colz");
555 
556  h_MuonTrackStitching_FailedCriteria_CriteriaTwoTracks_Angle = tfs->make<TH2D>("h_MuonTrackStitching_FailedCriteria_CriteriaTwoTracks_Angle", "Bad events: CriteriaTwoTracks vs. angle; Criteria; Angle [#circ]",7,0.75,4.25,180,0,180.);
557  h_MuonTrackStitching_FailedCriteria_CriteriaTwoTracks_Angle->SetOption("colz");
558 
559  //Stitching variables: good events
560  h_MuonTrackStitching_MatchedCriteria_TrackRes_Completeness = tfs->make<TH2D>("h_MuonTrackStitching_MatchedCriteria_TrackRes_Completeness", "Good events: L_{reco}/L_{truth} (leading) vs. Completeness (leading); L_{reco}/L_{truth} (leading); Completeness (leading)",150,0.,1.5,150,0.,1.5);
561  h_MuonTrackStitching_MatchedCriteria_TrackRes_Completeness->SetOption("colz");
562 
563  h_MuonTrackStitching_MatchedCriteria_Distance_Angle = tfs->make<TH2D>("h_MuonTrackStitching_MatchedCriteria_Distance_Angle", "Good events: distance vs. angle b/w leading and second muon track; Distance [cm]; Angle [#circ]",100,0.,100.,100,0.,180.);
564  h_MuonTrackStitching_MatchedCriteria_Distance_Angle->SetOption("colz");
565 
566  h_MuonTrackStitching_MatchedCriteria_TrackResLeading_TrackResSecond = tfs->make<TH2D>("h_MuonTrackStitching_MatchedCriteria_TrackResLeading_TrackResSecond", "Good events: L_{reco}/L_{truth} (leading) vs. L_{reco}/L_{truth} (second); L_{reco}/L_{truth} (leading); L_{reco}/L_{truth} (second)",150,0.,1.5,150,0.,1.5);
567  h_MuonTrackStitching_MatchedCriteria_TrackResLeading_TrackResSecond->SetOption("colz");
568 
569  h_MuonTrackStitching_MatchedCriteria_CriteriaTwoTracks_Angle = tfs->make<TH2D>("h_MuonTrackStitching_MatchedCriteria_CriteriaTwoTracks_Angle", "Good events: CriteriaTwoTracks vs. angle b/w leading and second muon track; Criteria; Angle [#circ]",7,0.75,4.25,100,0.,180.);
570  h_MuonTrackStitching_MatchedCriteria_CriteriaTwoTracks_Angle->SetOption("colz");
571 
572  //Stitching variables: bad events but leading plus second ok
573  h_MuonTrackStitching_FailedCriteriaButLeadingPlusSecondGood_TrackRes_Completeness = tfs->make<TH2D>("h_MuonTrackStitching_FailedCriteriaButLeadingPlusSecondGood_TrackRes_Completeness", "Bad events but leading + second is good: L_{reco}/L_{truth} (leading) vs. Completeness (leading); L_{reco}/L_{truth} (leading); Completeness (leading)",150,0.,1.5,150,0.,1.5);
574  h_MuonTrackStitching_FailedCriteriaButLeadingPlusSecondGood_TrackRes_Completeness->SetOption("colz");
575 
576  h_MuonTrackStitching_FailedCriteriaButLeadingPlusSecondGood_Distance_Angle = tfs->make<TH2D>("h_MuonTrackStitching_FailedCriteriaButLeadingPlusSecondGood_Distance_Angle", "Bad events but leading + second is good: distance vs. angle b/w leading and second muon track; Distance [cm]; Angle [#circ]",100,0.,100.,100,0.,180.);
577  h_MuonTrackStitching_FailedCriteriaButLeadingPlusSecondGood_Distance_Angle->SetOption("colz");
578 
579  h_MuonTrackStitching_FailedCriteriaButLeadingPlusSecondGood_TrackResLeading_TrackResSecond = tfs->make<TH2D>("h_MuonTrackStitching_FailedCriteriaButLeadingPlusSecondGood_TrackResLeading_TrackResSecond", "Bad events but leading + second is good: L_{reco}/L_{truth} (leading) vs. L_{reco}/L_{truth} (second); L_{reco}/L_{truth} (leading); L_{reco}/L_{truth} (second)",150,0.,1.5,150,0.,1.5);
580  h_MuonTrackStitching_FailedCriteriaButLeadingPlusSecondGood_TrackResLeading_TrackResSecond->SetOption("colz");
581 
582  //Stitching variables: bad events and leading plus second not ok
583  h_MuonTrackStitching_FailedCriteriaAndLeadingPlusSecondBad_TrackRes_Completeness = tfs->make<TH2D>("h_MuonTrackStitching_FailedCriteriaAndLeadingPlusSecondBad_TrackRes_Completeness", "Bad events and leading + second is bad: L_{reco}/L_{truth} (leading) vs. Completeness (leading); L_{reco}/L_{truth} (leading); Completeness (leading)",150,0.,1.5,150,0.,1.5);
584  h_MuonTrackStitching_FailedCriteriaAndLeadingPlusSecondBad_TrackRes_Completeness->SetOption("colz");
585 
586  h_MuonTrackStitching_FailedCriteriaAndLeadingPlusSecondBad_Distance_Angle = tfs->make<TH2D>("h_MuonTrackStitching_FailedCriteriaAndLeadingPlusSecondBad_Distance_Angle", "Bad events and leading + second is bad: distance vs. angle b/w leading and second muon track; Distance [cm]; Angle [#circ]",100,0.,100.,100,0.,180.);
587  h_MuonTrackStitching_FailedCriteriaAndLeadingPlusSecondBad_Distance_Angle->SetOption("colz");
588 
589  h_MuonTrackStitching_FailedCriteriaAndLeadingPlusSecondBad_TrackResLeading_TrackResSecond = tfs->make<TH2D>("h_MuonTrackStitching_FailedCriteriaAndLeadingPlusSecondBad_TrackResLeading_TrackResSecond", "Bad events and leading + second is bad: L_{reco}/L_{truth} (leading) vs. L_{reco}/L_{truth} (second); L_{reco}/L_{truth} (leading); L_{reco}/L_{truth} (second)",150,0.,1.5,150,0.,1.5);
590  h_MuonTrackStitching_FailedCriteriaAndLeadingPlusSecondBad_TrackResLeading_TrackResSecond->SetOption("colz");
591 
592 
593 
594 }
595 //========================================================================
596 void MuonTrackingEff::endJob(){
597 
598  doEfficiencies();
599 
600 }
601 //========================================================================
602 void MuonTrackingEff::beginRun(const art::Run& /*run*/){
603  mf::LogInfo("MuonTrackingEff")<<"begin run..."<<std::endl;
604 }
605 //========================================================================
606 void MuonTrackingEff::analyze( const art::Event& event ){
607  if (event.isRealData()) return;
608 
609  bool isFiducial = false;
610  processEff(event, isFiducial);
611 }
612 //========================================================================
613 void MuonTrackingEff::processEff( const art::Event& event, bool &isFiducial){
614 
615  EventCounter++;
616  simb::MCParticle *MCTruthMuonParticle = NULL;
617 
618  art::ServiceHandle<cheat::ParticleInventoryService> pi_serv;
619  const sim::ParticleList& plist = pi_serv->ParticleList();
620  simb::MCParticle *particle=0;
621 
622  for( sim::ParticleList::const_iterator ipar = plist.begin(); ipar!=plist.end(); ++ipar){
623  particle = ipar->second;
624 
625  if( particle->Mother() == 0 ){ //0=particle has no mother particle, 1=particle has a mother particle
626  const TLorentzVector& positionStart = particle->Position(0);
627  positionStart.GetXYZT(MCTruthMuonVertex); //MCTruthMuonVertex[0-2]=truth vertex, MCTruthMuonVertex[3]=t=0
628  }
629 
630  if( particle->PdgCode() == fMuonPDGCode ){ // Particle cannon muon
631  MCTruthMuonParticle = particle;
632  MCTruthMuonID = particle->TrackId();
633  MCTruthMuonMomentum = sqrt(pow(particle->Momentum().Px(),2)+pow(particle->Momentum().Py(),2)+pow(particle->Momentum().Pz(),2));
634 
635  if(particle->Momentum().Pz() >=0 && particle->Momentum().Px() >=0)
636  {
637  MCTruthMuonThetaXZ = (180.0/3.14159)*atan(particle->Momentum().Px()/particle->Momentum().Pz());
638  }
639  else if(particle->Momentum().Pz() < 0 && particle->Momentum().Px() >=0)
640  {
641  MCTruthMuonThetaXZ = 180.0 + (180.0/3.14159)*atan(particle->Momentum().Px()/particle->Momentum().Pz());
642  }
643  else if(particle->Momentum().Pz() < 0 && particle->Momentum().Px() < 0)
644  {
645  MCTruthMuonThetaXZ = 180.0 + (180.0/3.14159)*atan(particle->Momentum().Px()/particle->Momentum().Pz());
646  }
647  else if(particle->Momentum().Pz() >= 0 && particle->Momentum().Px() < 0)
648  {
649  MCTruthMuonThetaXZ = 360.0 + (180.0/3.14159)*atan(particle->Momentum().Px()/particle->Momentum().Pz());
650  }
651 
652  MCTruthMuonThetaYZ = (180.0/3.14159)*asin(particle->Momentum().Py()/MCTruthMuonMomentum);
653  }
654  }
655  double MCTruthLengthMuon = truthLength(MCTruthMuonParticle);
656  h_TruthLength->Fill(MCTruthLengthMuon);
657 
658  //===================================================================
659  //Saving denominator histograms
660  //===================================================================
661  isFiducial =insideFV( MCTruthMuonVertex );
662  if( !isFiducial ) return;
663 
664  //save events for Nucleon decay and particle cannon
665  if( MCTruthMuonParticle ){
666  h_ThetaXZ_den->Fill(MCTruthMuonThetaXZ);
667  h_ThetaYZ_den->Fill(MCTruthMuonThetaYZ);
668  h_SinThetaYZ_den->Fill(sin((3.14159/180.)*MCTruthMuonThetaYZ));
669  h_ThetaXZ_ThetaYZ_den->Fill(MCTruthMuonThetaXZ,MCTruthMuonThetaYZ);
670  h_ThetaXZ_SinThetaYZ_den->Fill(MCTruthMuonThetaXZ,sin((3.14159/180.)*MCTruthMuonThetaYZ));
671  CountMCTruthMuon++;
672  }
673 
674 
675 
676  //========================================================================
677  // Reco stuff, once we have selected a MC Particle let's find out if there is a track associated
678  //========================================================================
679 
680  int NMuonTracks=0;
681 
682  art::Handle<std::vector<recob::Track>> TrackListHandle;
683  if(! event.getByLabel(fTrackModuleLabel, TrackListHandle)) return;
684  std::vector<art::Ptr<recob::Track> > TrackList;
685  art::fill_ptr_vector(TrackList, TrackListHandle);
686  int NRecoTracks = TrackList.size();
687  art::FindManyP<recob::Hit> track_hits(TrackListHandle, event, fTrackModuleLabel);
688  if( NRecoTracks == 0 ){
689  LOG_DEBUG("MuonTrackingEff")<<"There are no reco tracks... bye";
690  std::cout << "There are no reco tracks! MCTruthMuonThetaXZ: " << std::endl;
691 
692  Criteria=0.;
693  h_FailedReconstruction_ThetaXZ_ThetaYZ->Fill(MCTruthMuonThetaXZ,MCTruthMuonThetaYZ);
694  h_FailedReconstruction_ThetaXZ_SinThetaYZ->Fill(MCTruthMuonThetaXZ,sin((3.14159/180.)*MCTruthMuonThetaYZ));
695  h_NoRecoTrackAtAll_ThetaXZ_SinThetaYZ->Fill(MCTruthMuonThetaXZ,sin((3.14159/180.)*MCTruthMuonThetaYZ));
696  h_Criteria_NRecoTrack_num->Fill(Criteria,static_cast<double>(NRecoTracks));
697  h_Criteria_NMuonTrack_num->Fill(Criteria,static_cast<double>(NMuonTracks));
698  CountNoRecoTracks++;
699  return;
700  }
701  LOG_DEBUG("MuonTrackingEff")<<"Found this many reco tracks "<<NRecoTracks;
702 
703  //std::cout << "NRecoTracks: " << NRecoTracks << std::endl;
704 
705  double PurityLeadingMuon=0.;
706  double CompletenessLeadingMuon=0.;
707  double RecoLengthLeadingMuon=0.;
708  art::Ptr<recob::Track> TrackLeadingMuon;
709 
710  double RecoLengthSecondMuon=0.;
711  double CompletenessSecondMuon=0.;
712  double PuritySecondMuon=0.;
713  art::Ptr<recob::Track> TrackSecondMuon;
714 
715  double TrackLengthMuonSum=0.;
716  double tmpTotalRecoEnergy=0.;
717 
718  double MaxLengthNoRecoMuon=0;
719  int PDGCodeMaxLengthNoRecoMuon=0;
720 
721  const simb::MCParticle *RecoMuonParticle = NULL;
722 
723  std::vector<art::Ptr<recob::Hit>> tmp_TrackHits = track_hits.at(0);
724  std::vector<art::Ptr<recob::Hit>> AllHits;
725  art::Handle<std::vector<recob::Hit>> HitHandle;
726  if(event.get(tmp_TrackHits[0].id(), HitHandle)) art::fill_ptr_vector(AllHits, HitHandle);
727 
728 
729  //Loop over reco tracks
730  for(int i=0; i<NRecoTracks; i++) {
731  art::Ptr<recob::Track> track = TrackList[i];
732  std::vector<art::Ptr<recob::Hit>> TrackHits = track_hits.at(i);
733  double tmpPurity=0.;
734  double tmpCompleteness=0.;
735  const simb::MCParticle *particle;
736 
737  truthMatcher(AllHits, TrackHits, particle, tmpPurity, tmpCompleteness, tmpTotalRecoEnergy);
738 
739  if (!particle) { std::cout << "ERROR: Truth matcher didn't find a particle!" << std::endl; continue;}
740 
741  if(track->Length() > MaxLengthNoRecoMuon && particle->PdgCode() != fMuonPDGCode && particle->TrackId() != MCTruthMuonID)
742  {
743  MaxLengthNoRecoMuon = track->Length();
744  PDGCodeMaxLengthNoRecoMuon = particle->PdgCode();
745  }
746  //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
747  if( (particle->PdgCode() == fMuonPDGCode) && (particle->TrackId() == MCTruthMuonID) && tmpCompleteness > 0 && tmpPurity > 0 ){
748 
749  NMuonTracks++;
750  TrackLengthMuonSum+=track->Length();
751 
752  if(NMuonTracks==1)
753  {
754  CompletenessLeadingMuon = tmpCompleteness;
755  PurityLeadingMuon = tmpPurity;
756  RecoLengthLeadingMuon = track->Length();
757  TrackLeadingMuon = track;
758 
759  RecoMuonParticle = particle;
760  }
761 
762 
763  if(NMuonTracks>=2 && tmpCompleteness > CompletenessLeadingMuon )
764  {
765 
766  CompletenessSecondMuon = CompletenessLeadingMuon;
767  PuritySecondMuon = PurityLeadingMuon;
768  RecoLengthSecondMuon = RecoLengthLeadingMuon;
769  TrackSecondMuon = TrackLeadingMuon;
770 
771  CompletenessLeadingMuon = tmpCompleteness;
772  PurityLeadingMuon = tmpPurity;
773  RecoLengthLeadingMuon = track->Length();
774  TrackLeadingMuon = track;
775 
776  RecoMuonParticle = particle;
777  }
778 
779  else if(NMuonTracks>=2 && tmpCompleteness < CompletenessLeadingMuon && tmpCompleteness > CompletenessSecondMuon)
780  {
781  CompletenessSecondMuon = tmpCompleteness;
782  PuritySecondMuon = tmpPurity;
783  RecoLengthSecondMuon = track->Length();
784  TrackSecondMuon = track;
785  }
786  }
787  }
788 
789  h_TotalRecoEnergy->Fill(tmpTotalRecoEnergy);
790 
791 
792 
793  //Muon events
794  //=========================================================
795 
796  if( RecoMuonParticle && MCTruthMuonParticle )
797  {
798  CountRecoMuon++;
799  h_Purity->Fill(PurityLeadingMuon);
800  h_Completeness->Fill(CompletenessLeadingMuon);
801  h_TrackRes->Fill(RecoLengthLeadingMuon/MCTruthLengthMuon);
802 
803  std::cout << "TrackLeadingMuon->Vertex().X(): " << TrackLeadingMuon->Vertex().X() << std::endl;
804  std::cout << "MCTruthMuonParticle->Vz(): " << MCTruthMuonParticle->Vz() << std::endl;
805  std::cout << " " << std::endl;
806 
807  double DistanceBetweenTruthAndRecoTrack;
808  double AngleBetweenTruthAndRecoTrack;
809  FuncDistanceAndAngleBetweenTruthAndRecoTrack(RecoMuonParticle,TrackLeadingMuon,DistanceBetweenTruthAndRecoTrack,AngleBetweenTruthAndRecoTrack);
810 
811  h_VertexRes->Fill(DistanceBetweenTruthAndRecoTrack);
812  h_DirectionRes->Fill(AngleBetweenTruthAndRecoTrack);
813 
814  h_MuonTrackStitching_TrackRes_Completeness->Fill(RecoLengthLeadingMuon/MCTruthLengthMuon, CompletenessLeadingMuon);
815  if(NMuonTracks>=2)
816  {
817  double DistanceBetweenTracks;
818  double AngleBetweenTracks;
819  double CriteriaTwoTracks;
820 
821  FuncDistanceAndAngleBetweenTracks(TrackSecondMuon,TrackLeadingMuon,DistanceBetweenTracks, AngleBetweenTracks, CriteriaTwoTracks);
822 
823  h_MuonTrackStitching_TrackResLeading_TrackResSecond->Fill(RecoLengthLeadingMuon/MCTruthLengthMuon,RecoLengthSecondMuon/MCTruthLengthMuon);
824  h_MuonTrackStitching_Distance_Angle->Fill(DistanceBetweenTracks,AngleBetweenTracks);
825  h_MuonTrackStitching_TrackResSecondMuon_Angle->Fill(RecoLengthSecondMuon/MCTruthLengthMuon,AngleBetweenTracks);
826  h_MuonTrackStitching_CompletenessSecondMuon_Angle->Fill(CompletenessSecondMuon,AngleBetweenTracks);
827  h_MuonTrackStitching_CriteriaTwoTracks_Angle->Fill(CriteriaTwoTracks, AngleBetweenTracks);
828  }
829 
830  //Completeness
831  if(CompletenessLeadingMuon < 0.5)
832  {
833  CountCompleteness++;
834  Criteria=2.;
835  h_Completeness_ThetaXZ_SinThetaYZ->Fill(MCTruthMuonThetaXZ,sin((3.14159/180.)*MCTruthMuonThetaYZ));
836  h_Criteria_NRecoTrack_num->Fill(Criteria,static_cast<double>(NRecoTracks));
837  h_Criteria_NMuonTrack_num->Fill(Criteria,static_cast<double>(NMuonTracks));
838  }
839  //Purity
840  if(PurityLeadingMuon < 0.5)
841  {
842  CountPurity++;
843  Criteria=3.;
844  h_Purity_ThetaXZ_SinThetaYZ->Fill(MCTruthMuonThetaXZ,sin((3.14159/180.)*MCTruthMuonThetaYZ));
845  h_Criteria_NRecoTrack_num->Fill(Criteria,static_cast<double>(NRecoTracks));
846  h_Criteria_NMuonTrack_num->Fill(Criteria,static_cast<double>(NMuonTracks));
847  }
848  //Track too short
849  if(RecoLengthLeadingMuon/MCTruthLengthMuon < 0.75)
850  {
851  CountTrackLengthTooShort++;
852  Criteria=4.;
853  h_TrackTooShort_ThetaXZ_SinThetaYZ->Fill(MCTruthMuonThetaXZ,sin((3.14159/180.)*MCTruthMuonThetaYZ));
854  h_Criteria_NRecoTrack_num->Fill(Criteria,static_cast<double>(NRecoTracks));
855  h_Criteria_NMuonTrack_num->Fill(Criteria,static_cast<double>(NMuonTracks));
856  }
857  //Track too long
858  if(RecoLengthLeadingMuon/MCTruthLengthMuon > 1.25)
859  {
860  CountTrackLengthTooLong++;
861  Criteria=5.;
862  h_TrackTooLong_ThetaXZ_SinThetaYZ->Fill(MCTruthMuonThetaXZ,sin((3.14159/180.)*MCTruthMuonThetaYZ));
863  h_Criteria_NRecoTrack_num->Fill(Criteria,static_cast<double>(NRecoTracks));
864  h_Criteria_NMuonTrack_num->Fill(Criteria,static_cast<double>(NMuonTracks));
865 
866  /*std::cout << "Track too long!" << std::endl;
867  art::ServiceHandle<cheat::BackTracker> bt2;
868  const sim::ParticleList& plist2 = bt2->ParticleList();
869  simb::MCParticle *particle2=0;
870 
871  std::cout << "EventCounter: " << EventCounter << std::endl;
872  for( sim::ParticleList::const_iterator ipar2 = plist2.begin(); ipar2!=plist2.end(); ++ipar2)
873  {
874  particle2 = ipar2->second;
875  std::cout << "particle2->TrackId(): " << particle2->TrackId() << std::endl;
876  std::cout << "particle2->PdgCode(): " << particle2->PdgCode() << std::endl;
877  std::cout << "particle2->Momentum().P(): " << particle2->Momentum().P() << std::endl;
878  std::cout << "tuthLength(particle2): " << truthLength(particle2) << std::endl;
879  std::cout << "particle2->Position(): (x,y,z): " << particle2->Vx() << "\t" << particle2->Vy() << "\t" << particle2->Vz() << std::endl;
880  std::cout << "particle2->Momentum(): (Px,Py,Pz): " << particle2->Momentum().Px() << "\t" << particle2->Momentum().Py() << "\t" << particle2->Momentum().Pz() << std::endl;
881  std::cout << "particle2->Position().T(): " << particle2->Position().T() << std::endl;
882  std::cout << "" << std::endl;
883  }*/
884 
885  }
886  //Reco failed at least one of the above criteria
887  if(CompletenessLeadingMuon < 0.5 || PurityLeadingMuon < 0.5 || RecoLengthLeadingMuon/MCTruthLengthMuon < 0.75 || RecoLengthLeadingMuon/MCTruthLengthMuon > 1.25)
888  {
889  CountBadLeadingMuonTrack++;
890  h_FailedReconstruction_ThetaXZ_ThetaYZ->Fill(MCTruthMuonThetaXZ,MCTruthMuonThetaYZ);
891  h_FailedReconstruction_ThetaXZ_SinThetaYZ->Fill(MCTruthMuonThetaXZ,sin((3.14159/180.)*MCTruthMuonThetaYZ));
892  h_MuonTrackStitching_FailedCriteria_TrackRes_Completeness->Fill(RecoLengthLeadingMuon/MCTruthLengthMuon, CompletenessLeadingMuon);
893 
894 
895  if(NMuonTracks==1) {BadEvents1MuonTrack++;}
896  if(NMuonTracks==2) {BadEvents2MuonTrack++;}
897  if(NMuonTracks==3) {BadEvents3MuonTrack++;}
898  if(NMuonTracks==4) {BadEvents4OrMoreMuonTrack++;}
899 
900  if(NMuonTracks>=2)
901  {
902  double AngleBetweenTracks;
903  double DistanceBetweenTracks;
904  double CriteriaTwoTracks;
905  FuncDistanceAndAngleBetweenTracks(TrackSecondMuon,TrackLeadingMuon,DistanceBetweenTracks, AngleBetweenTracks, CriteriaTwoTracks);
906 
907  if(AngleBetweenTracks>160.)
908  {
909  std::cout << "EventCounter: " << EventCounter << std::endl;
910  std::cout << "Angle b/w tracks: " << AngleBetweenTracks << std::endl;
911  std::cout << "CriteriaTwoTracks: " << CriteriaTwoTracks << std::endl;
912  std::cout << "CompletenessLeadingMuon: " << CompletenessLeadingMuon << std::endl;
913  std::cout << "CompletenessSecondMuon: " << CompletenessSecondMuon << std::endl;
914  std::cout << "PurityLeadingMuon: " << PurityLeadingMuon << std::endl;
915  std::cout << "PuritySecondMuon: " << PuritySecondMuon << std::endl;
916  std::cout << "TrackLeadingMuon: " << RecoLengthLeadingMuon/MCTruthLengthMuon << std::endl;
917  std::cout << "TrackResSecondMuon: " << RecoLengthSecondMuon/MCTruthLengthMuon << std::endl;
918  std::cout << "TrackID leading: " << TrackLeadingMuon->ID() << std::endl;
919  std::cout << "TrackID second: " << TrackSecondMuon->ID() << std::endl;
920  }
921 
922  h_MuonTrackStitching_FailedCriteria_Distance_Angle->Fill(DistanceBetweenTracks,AngleBetweenTracks);
923  h_MuonTrackStitching_FailedCriteria_TrackResLeading_TrackResSecond->Fill(RecoLengthLeadingMuon/MCTruthLengthMuon,RecoLengthSecondMuon/MCTruthLengthMuon);
924  h_MuonTrackStitching_FailedCriteria_CompletenessLeading_CompletenessSecond->Fill(CompletenessLeadingMuon,CompletenessSecondMuon);
925  h_MuonTrackStitching_FailedCriteria_TrackResSecondMuon_Angle->Fill(RecoLengthSecondMuon/MCTruthLengthMuon, AngleBetweenTracks);
926  h_MuonTrackStitching_FailedCriteria_CompletenessSecondMuon_Angle->Fill(CompletenessSecondMuon, AngleBetweenTracks);
927  h_MuonTrackStitching_FailedCriteria_CriteriaTwoTracks_Angle->Fill(CriteriaTwoTracks, AngleBetweenTracks);
928  if( (CompletenessLeadingMuon+CompletenessSecondMuon) >= 0.5 && PurityLeadingMuon >= 0.5 && PuritySecondMuon >= 0.5 && (RecoLengthLeadingMuon+RecoLengthSecondMuon)/MCTruthLengthMuon >= 0.75 && (RecoLengthLeadingMuon+RecoLengthSecondMuon)/MCTruthLengthMuon <= 1.25 ) {
929  CountBadLeadingMuonTrackButLeadingPlusSecondGood++;
930  h_ThetaXZ_ThetaYZ_LeadingPlusSecondOk->Fill(MCTruthMuonThetaXZ,MCTruthMuonThetaYZ);
931  h_ThetaXZ_SinThetaYZ_LeadingPlusSecondOk->Fill(MCTruthMuonThetaXZ,sin((3.14159/180.)*MCTruthMuonThetaYZ));
932 
933  h_MuonTrackStitching_FailedCriteriaButLeadingPlusSecondGood_TrackRes_Completeness->Fill(RecoLengthLeadingMuon/MCTruthLengthMuon, CompletenessLeadingMuon);
934  h_MuonTrackStitching_FailedCriteriaButLeadingPlusSecondGood_TrackResLeading_TrackResSecond->Fill(RecoLengthLeadingMuon/MCTruthLengthMuon,RecoLengthSecondMuon/MCTruthLengthMuon);
935  h_MuonTrackStitching_FailedCriteriaButLeadingPlusSecondGood_Distance_Angle->Fill(DistanceBetweenTracks,AngleBetweenTracks);
936  }
937  if (( CompletenessLeadingMuon+CompletenessSecondMuon) < 0.5 || PurityLeadingMuon < 0.5 || PuritySecondMuon < 0.5 || (RecoLengthLeadingMuon+RecoLengthSecondMuon)/MCTruthLengthMuon < 0.75 || (RecoLengthLeadingMuon+RecoLengthSecondMuon)/MCTruthLengthMuon > 1.25 )
938  {
939  CountBadLeadingMuonTrackAndLeadingPlusSecondBad++;
940 
941  h_MuonTrackStitching_FailedCriteriaAndLeadingPlusSecondBad_TrackRes_Completeness->Fill(RecoLengthLeadingMuon/MCTruthLengthMuon, CompletenessLeadingMuon);
942  h_MuonTrackStitching_FailedCriteriaAndLeadingPlusSecondBad_TrackResLeading_TrackResSecond->Fill(RecoLengthLeadingMuon/MCTruthLengthMuon,RecoLengthSecondMuon/MCTruthLengthMuon);
943  h_MuonTrackStitching_FailedCriteriaAndLeadingPlusSecondBad_Distance_Angle->Fill(DistanceBetweenTracks,AngleBetweenTracks);
944  }
945  if((CompletenessLeadingMuon+CompletenessSecondMuon) < 0.5) {CountBadLeadingMuonTrackAndLeadingPlusSecondBadCompleteness++;}
946  if(PurityLeadingMuon < 0.5 || PuritySecondMuon < 0.5 ) {CountBadLeadingMuonTrackAndLeadingPlusSecondBadPurity++;}
947  if((RecoLengthLeadingMuon+RecoLengthSecondMuon)/MCTruthLengthMuon < 0.75) {CountBadLeadingMuonTrackAndLeadingPlusSecondBadTrackTooShort++;}
948  if((RecoLengthLeadingMuon+RecoLengthSecondMuon)/MCTruthLengthMuon > 1.25) {CountBadLeadingMuonTrackAndLeadingPlusSecondBadTrackTooLong++;}
949 
950  }
951  else if(NMuonTracks==1)
952  {
953  CountBadLeadingMuonTrackAndOnlyOneMuonTrack++;
954  if(CompletenessLeadingMuon < 0.5) {CountBadLeadingMuonTrackAndOnlyOneMuonTrackCompleteness++;}
955  if(PurityLeadingMuon < 0.5) {CountBadLeadingMuonTrackAndOnlyOneMuonTrackPurity++;}
956  if(RecoLengthLeadingMuon/MCTruthLengthMuon < 0.75) {CountBadLeadingMuonTrackAndOnlyOneMuonTrackTrackTooShort++;}
957  if(RecoLengthLeadingMuon/MCTruthLengthMuon > 1.25) {CountBadLeadingMuonTrackAndOnlyOneMuonTrackTrackTooLong++;}
958  }
959  }
960  //Reco succesful according to the above criteria
961  if(CompletenessLeadingMuon >= 0.5 && PurityLeadingMuon >= 0.5 && RecoLengthLeadingMuon/MCTruthLengthMuon >= 0.75 && RecoLengthLeadingMuon/MCTruthLengthMuon <= 1.25)
962  {
963  CountGoodLeadingMuonTrack++;
964  Criteria=6.;
965  h_ThetaXZ_num->Fill(MCTruthMuonThetaXZ);
966  h_ThetaYZ_num->Fill(MCTruthMuonThetaYZ);
967  h_SinThetaYZ_num->Fill(sin((3.14159/180.)*MCTruthMuonThetaYZ));
968  h_ThetaXZ_ThetaYZ_num->Fill(MCTruthMuonThetaXZ,MCTruthMuonThetaYZ);
969  h_ThetaXZ_SinThetaYZ_num->Fill(MCTruthMuonThetaXZ,sin((3.14159/180.)*MCTruthMuonThetaYZ));
970  h_Criteria_NRecoTrack_num->Fill(Criteria,static_cast<double>(NRecoTracks));
971  h_Criteria_NMuonTrack_num->Fill(Criteria,static_cast<double>(NMuonTracks));
972 
973  h_MuonTrackStitching_MatchedCriteria_TrackRes_Completeness->Fill(RecoLengthLeadingMuon/MCTruthLengthMuon, CompletenessLeadingMuon);
974 
975  if(NMuonTracks==1) {GoodEvents1MuonTrack++;}
976  if(NMuonTracks==2) {GoodEvents2MuonTrack++;}
977  if(NMuonTracks==3) {GoodEvents3MuonTrack++;}
978  if(NMuonTracks==4) {GoodEvents4OrMoreMuonTrack++;}
979 
980  if(NMuonTracks>=2)
981  {
982  h_MuonTrackStitching_MatchedCriteria_TrackResLeading_TrackResSecond->Fill(RecoLengthLeadingMuon/MCTruthLengthMuon,RecoLengthSecondMuon/MCTruthLengthMuon);
983 
984  double AngleBetweenTracks;
985  double DistanceBetweenTracks;
986  double CriteriaTwoTracks;
987  FuncDistanceAndAngleBetweenTracks(TrackSecondMuon,TrackLeadingMuon,DistanceBetweenTracks, AngleBetweenTracks, CriteriaTwoTracks);
988  h_MuonTrackStitching_MatchedCriteria_Distance_Angle->Fill(DistanceBetweenTracks,AngleBetweenTracks);
989  h_MuonTrackStitching_MatchedCriteria_CriteriaTwoTracks_Angle->Fill(CriteriaTwoTracks, AngleBetweenTracks);
990  }
991  }
992  }
993  //No muon track
994  if( !RecoMuonParticle && MCTruthMuonParticle )
995  {
996  CountNoMuonTracks++;
997  BadEvents0MuonTrack++;
998  Criteria=1.;
999  h_NoMuonTrack_ThetaXZ_SinThetaYZ->Fill(MCTruthMuonThetaXZ,sin((3.14159/180.)*MCTruthMuonThetaYZ));
1000  h_FailedReconstruction_ThetaXZ_ThetaYZ->Fill(MCTruthMuonThetaXZ,MCTruthMuonThetaYZ);
1001  h_FailedReconstruction_ThetaXZ_SinThetaYZ->Fill(MCTruthMuonThetaXZ,sin((3.14159/180.)*MCTruthMuonThetaYZ));
1002  h_Criteria_NRecoTrack_num->Fill(Criteria,static_cast<double>(NRecoTracks));
1003  h_Criteria_NMuonTrack_num->Fill(Criteria,static_cast<double>(NMuonTracks));
1004  h_NoMuonTrack_MaxTrackLength_PDGCode->Fill(MaxLengthNoRecoMuon, static_cast<double>(PDGCodeMaxLengthNoRecoMuon));
1005  }
1006 }
1007 //========================================================================
1008 void MuonTrackingEff::truthMatcher( 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){
1009 
1010  //std::cout<<"truthMatcher..."<<std::endl;
1011  art::ServiceHandle<cheat::BackTrackerService> bt_serv;
1012  art::ServiceHandle<cheat::ParticleInventoryService> pi_serv;
1013  std::map<int,double> trkID_E; //map that connects TrackID and energy for each hit <trackID, energy>
1014  for(size_t j = 0; j < track_hits.size(); ++j){ //loop over all the hits in this track
1015  art::Ptr<recob::Hit> hit = track_hits[j];
1016  std::vector<sim::TrackIDE> TrackIDs = bt_serv->HitToTrackIDEs(hit); //TrackIDE contains TrackID, energy and energyFrac. A hit can have several TrackIDs (so this hit is associated with multiple MC truth track IDs (EM shower IDs are negative). If a hit ahs multiple trackIDs, "energyFrac" contains the fraction of the energy of for each ID compared to the total energy of the hit. "energy" contains only the energy associated with the specific ID in that case. This requires MC truth info!
1017  for(size_t k = 0; k < TrackIDs.size(); k++){ //Loop over the TrackIDs of each hit
1018  trkID_E[TrackIDs[k].trackID] += TrackIDs[k].energy; //sum up the energy for each TrackID and store <TrackID, energy> in "TrkID_E"
1019  }
1020  }
1021 
1022 
1023  double E_em =0.0;
1024  double max_E = -999.0;
1025  double TotalEnergyTrack = 0.0;
1026  int TrackID = -999;
1027  double PartialEnergyTrackID=0.0; // amount of energy deposited by the particle that deposited more energy... tomato potato... blabla
1030  if( !trkID_E.size() ) {
1031  MCparticle = 0;
1032  return; //Ghost track???
1033  }
1034  for(std::map<int,double>::iterator ii = trkID_E.begin(); ii!=trkID_E.end(); ++ii){ //trkID_E contains the trackID (first) and corresponding energy (second) for a specific track, summed up over all events. here looping over all trekID_E's
1035  TotalEnergyTrack += ii->second; //and summing up the energy of all hits in the track (TotalEnergyTrack)
1036  if((ii->second)>max_E){ //looking for the trakID with the highest energy in the track. this is PartialEnergyTrackID and max_E then
1037  PartialEnergyTrackID = ii->second;
1038  max_E = ii->second;
1039  TrackID = ii->first; //saving trackID of the ID with the highest energy contribution in the track to assign it to MCparticle later
1040  if( TrackID < 0 ) E_em += ii->second; //IDs of em shower particles are negative
1041  }
1042  }
1043  MCparticle = pi_serv->TrackIdToParticle_P(TrackID);
1044  //In the current simulation, we do not save EM Shower daughters in GEANT. But we do save the energy deposition in TrackIDEs. If the energy deposition is from a particle that is the daughter of
1045  //an EM particle, the negative of the parent track ID is saved in TrackIDE for the daughter particle
1046  //we don't want to track gammas or any other EM activity
1047  if( TrackID < 0 ) {return;}
1048 
1049  //Purity = (PartialEnergyTrackID+E_em)/TotalEnergyTrack;
1050  Purity = PartialEnergyTrackID/TotalEnergyTrack;
1051 
1052  //completeness
1053  TotalRecoEnergy =0;
1054  for(size_t k = 0; k < AllHits.size(); ++k){ //loop over all hits (all hits in all tracks of the event, not only the hits in the track we were looking at before)
1055  art::Ptr<recob::Hit> hit = AllHits[k];
1056  std::vector<sim::TrackIDE> TrackIDs = bt_serv->HitToTrackIDEs(hit);
1057  for(size_t l = 0; l < TrackIDs.size(); ++l){ //and over all track IDs of the hits
1058  if(TrackIDs[l].trackID==TrackID) TotalRecoEnergy += TrackIDs[l].energy; //and sum up the energy fraction of all hits that correspond ot the saved trackID
1059  }
1060  }
1061  Completeness = PartialEnergyTrackID/TotalRecoEnergy;
1062 }
1063 
1064 void MuonTrackingEff::FuncDistanceAndAngleBetweenTracks(art::Ptr<recob::Track> Track1, art::Ptr<recob::Track> Track2, double &TempDistanceBetweenTracks, double &TempAngleBetweenTracks, double &TempCriteriaTwoTracks)
1065 {
1066 
1067 TempDistanceBetweenTracks = sqrt(pow(Track1->End().X()-Track2->Vertex().X(),2) + pow(Track1->End().Y()-Track2->Vertex().Y(),2) + pow(Track1->End().Z()-Track2->Vertex().Z(),2));
1068 TempAngleBetweenTracks = (180.0/3.14159)*Track1->EndDirection().Angle(Track2->VertexDirection());
1069 TempCriteriaTwoTracks = 1.;
1070 
1071 if(TempDistanceBetweenTracks > sqrt(pow(Track1->End().X()-Track2->End().X(),2) + pow(Track1->End().Y()-Track2->End().Y(),2) + pow(Track1->End().Z()-Track2->End().Z(),2)))
1072  {
1073  TempDistanceBetweenTracks = sqrt(pow(Track1->End().X()-Track2->End().X(),2) + pow(Track1->End().Y()-Track2->End().Y(),2) + pow(Track1->End().Z()-Track2->End().Z(),2));
1074  TempAngleBetweenTracks = 180. - (180.0/3.14159)*Track1->EndDirection().Angle(Track2->EndDirection());
1075  TempCriteriaTwoTracks = 2.;
1076  }
1077 
1078 if(TempDistanceBetweenTracks > sqrt(pow(Track1->Vertex().X()-Track2->End().X(),2) + pow(Track1->Vertex().Y()-Track2->End().Y(),2) + pow(Track1->Vertex().Z()-Track2->End().Z(),2)))
1079  {
1080  TempDistanceBetweenTracks = sqrt(pow(Track1->Vertex().X()-Track2->End().X(),2) + pow(Track1->Vertex().Y()-Track2->End().Y(),2) + pow(Track1->Vertex().Z()-Track2->End().Z(),2));
1081  TempAngleBetweenTracks = (180.0/3.14159)*Track1->VertexDirection().Angle(Track2->EndDirection());
1082  TempCriteriaTwoTracks = 3.;
1083  }
1084 
1085 if(TempDistanceBetweenTracks > sqrt(pow(Track1->Vertex().X()-Track2->Vertex().X(),2) + pow(Track1->Vertex().Y()-Track2->Vertex().Y(),2) + pow(Track1->Vertex().Z()-Track2->Vertex().Z(),2)))
1086  {
1087  TempDistanceBetweenTracks = sqrt(pow(Track1->Vertex().X()-Track2->Vertex().X(),2) + pow(Track1->Vertex().Y()-Track2->Vertex().Y(),2) + pow(Track1->Vertex().Z()-Track2->Vertex().Z(),2));
1088  TempAngleBetweenTracks = 180. - (180.0/3.14159)*Track1->VertexDirection().Angle(Track2->VertexDirection());
1089  TempCriteriaTwoTracks = 4.;
1090  }
1091 
1092 }
1093 
1094 void MuonTrackingEff::FuncDistanceAndAngleBetweenTruthAndRecoTrack(const simb::MCParticle *&MCparticle, art::Ptr<recob::Track> Track, double &TempDistanceBetweenTruthAndRecoTrack, double &TempAngleBeetweenTruthAndRecoTrack)
1095 {
1096 TempDistanceBetweenTruthAndRecoTrack = sqrt(pow(Track->Vertex().X() - MCparticle->Vx(),2) + pow(Track->Vertex().Y() - MCparticle->Vy(),2) + pow(Track->Vertex().Z() - MCparticle->Vz(),2));
1097 
1098 TempAngleBeetweenTruthAndRecoTrack=0;
1099 
1100 }
1101 
1102 //========================================================================
1103 double MuonTrackingEff::truthLength( const simb::MCParticle *MCparticle ){
1104  //calculate the truth length considering only the part that is inside the TPC
1105  //Base on a piece of code from dune/TrackingAna/TrackingEfficiency_module.cc
1106 
1107  if( !MCparticle ) return -999.0;
1108  int numberTrajectoryPoints = MCparticle->NumberTrajectoryPoints();
1109  std::vector<double> TPCLengthHits(numberTrajectoryPoints, 0);
1110  int FirstHit=0, LastHit=0;
1111  double TPCLength = 0.0;
1112  bool BeenInVolume = false;
1113 
1114  for(unsigned int MCHit=0; MCHit < TPCLengthHits.size(); ++MCHit) {
1115  const TLorentzVector& tmpPosition= MCparticle->Position(MCHit);
1116  double const tmpPosArray[]={tmpPosition[0],tmpPosition[1],tmpPosition[2]};
1117  if (MCHit!=0) TPCLengthHits[MCHit] = sqrt( pow( (MCparticle->Vx(MCHit-1)-MCparticle->Vx(MCHit)),2)+ pow( (MCparticle->Vy(MCHit-1)-MCparticle->Vy(MCHit)),2)+ pow( (MCparticle->Vz(MCHit-1)-MCparticle->Vz(MCHit)),2));
1118  geo::TPCID tpcid = geom->FindTPCAtPosition(tmpPosArray);
1119  if(tpcid.isValid) {
1120  // -- Check if hit is within drift window...
1121  /* geo::CryostatGeo const& cryo = geom->Cryostat(tpcid.Cryostat);
1122  geo::TPCGeo const& tpc = cryo.TPC(tpcid.TPC);
1123  double XPlanePosition = tpc.PlaneLocation(0)[0];
1124  std::cout << "XPlanePosition: " << XPlanePosition << std::endl;
1125 
1126  double DriftTimeCorrection = fabs( tmpPosition[0] - XPlanePosition ) / XDriftVelocity;
1127  std::cout << "DriftTimeCorrection: " << DriftTimeCorrection << std::endl;
1128  std::cout << "tmpPosition[0]: " << tmpPosition[0] << std::endl;
1129  std::cout << "XDriftVelocity: " << XDriftVelocity << std::endl;
1130  std::cout << "MCparticle->T(): " << MCparticle->T() << std::endl;
1131  double TimeAtPlane = MCparticle->T() + DriftTimeCorrection;
1132  std::cout << "TimeAtPlane: " << TimeAtPlane << "\t" << "detprop->TriggerOffset(): " << detprop->TriggerOffset() << std::endl;
1133  std::cout << "WindowSize: " << WindowSize << std::endl;
1134  if( TimeAtPlane < detprop->TriggerOffset() || TimeAtPlane > detprop->TriggerOffset() + WindowSize ){ std::cout << "BYE" << std::endl; continue;} */
1135  LastHit = MCHit;
1136  if( !BeenInVolume ) {
1137  BeenInVolume = true;
1138  FirstHit = MCHit;
1139  }
1140  }
1141  }
1142  for (int Hit = FirstHit+1; Hit <= LastHit; ++Hit ) TPCLength += TPCLengthHits[Hit];
1143  return TPCLength;
1144 }
1145 //========================================================================
1146 bool MuonTrackingEff::insideFV( double vertex[4]){
1147 
1148  double x = vertex[0];
1149  double y = vertex[1];
1150  double z = vertex[2];
1151 
1152  if (x>fFidVolXmin && x<fFidVolXmax&&
1153  y>fFidVolYmin && y<fFidVolYmax&&
1154  z>fFidVolZmin && z<fFidVolZmax)
1155  return true;
1156  else
1157  return false;
1158 }
1159 //========================================================================
1160 void MuonTrackingEff::doEfficiencies(){
1161  std::cout << std::endl;
1162 
1163  std::cout << "EventCounter: " << "\t" << EventCounter << std::endl;
1164 
1165  std::cout << "CountMCTruthMuon: " << "\t" << CountMCTruthMuon << " = " << 100*static_cast<double>(CountMCTruthMuon)/static_cast<double>(EventCounter) <<"%" << std::endl;
1166 
1167  std::cout << "CountGoodLeadingMuonTrack (=good events): " << "\t" << CountGoodLeadingMuonTrack << "/" << CountMCTruthMuon << " = " << 100*static_cast<double>(CountGoodLeadingMuonTrack)/static_cast<double>(CountMCTruthMuon) <<"%" << std::endl;
1168 
1169  std::cout << "CountBadLeadingMuonTrack+CountNoRecoTracks+CountNoMuonTracks (=bad events): " << "\t" << CountBadLeadingMuonTrack+CountNoRecoTracks+CountNoMuonTracks << " = " << 100*static_cast<double>(CountBadLeadingMuonTrack+CountNoRecoTracks+CountNoMuonTracks)/static_cast<double>(CountMCTruthMuon) <<"%" << std::endl;
1170 
1171  std::cout << "CountNoRecoTracks+CountNoMuonTracks: " << "\t" << CountNoRecoTracks+CountNoMuonTracks << " = " << 100*static_cast<double>(CountNoRecoTracks+CountNoMuonTracks)/static_cast<double>(CountMCTruthMuon) <<"%" << std::endl;
1172 
1173  std::cout << "CountTrackLengthTooShort: " << "\t" << CountTrackLengthTooShort << " = " << 100*static_cast<double>(CountTrackLengthTooShort)/static_cast<double>(CountMCTruthMuon) <<"%" << std::endl;
1174 
1175 
1176  std::cout << "CountCompleteness: " << "\t" << CountCompleteness << " = " << 100*static_cast<double>(CountCompleteness)/static_cast<double>(CountMCTruthMuon) <<"%" << std::endl;
1177 
1178  std::cout << "CountTrackLengthTooLong: " << "\t" << CountTrackLengthTooLong << " = " << 100*static_cast<double>(CountTrackLengthTooLong)/static_cast<double>(CountMCTruthMuon) <<"%" << std::endl;
1179 
1180  std::cout << "CountPurity: " << "\t" << CountPurity << " = " << 100*static_cast<double>(CountPurity)/static_cast<double>(CountMCTruthMuon) <<"%" << std::endl;
1181 
1182 
1183  std::cout << std::endl;
1184 
1185  std::cout << "GoodLeadingMuonTrack+CountBadLeadingMuonTrackButLeadingPlusSecondGood (=good events after stitching): " << "\t" << CountGoodLeadingMuonTrack+CountBadLeadingMuonTrackButLeadingPlusSecondGood << "/" << CountMCTruthMuon << " = " << 100*static_cast<double>(CountGoodLeadingMuonTrack+CountBadLeadingMuonTrackButLeadingPlusSecondGood)/static_cast<double>(CountMCTruthMuon) <<"%" << std::endl;
1186 
1187  std::cout << "CountBadLeadingMuonTrack+CountNoRecoTracks+CountNoMuonTracks-CountBadLeadingMuonTrackButLeadingPlusSecondGood (=bad events after stitching) : " << "\t" << CountBadLeadingMuonTrack+CountNoRecoTracks+CountNoMuonTracks-CountBadLeadingMuonTrackButLeadingPlusSecondGood << " = " << 100*static_cast<double>(CountBadLeadingMuonTrack+CountNoRecoTracks+CountNoMuonTracks-CountBadLeadingMuonTrackButLeadingPlusSecondGood)/static_cast<double>(CountMCTruthMuon) <<"%" << std::endl;
1188 
1189  std::cout << std::endl;
1190 
1191  std::cout << "CountBadLeadingMuonTrackButLeadingPlusSecondGood: " << "\t" << CountBadLeadingMuonTrackButLeadingPlusSecondGood << " = " << 100*static_cast<double>(CountBadLeadingMuonTrackButLeadingPlusSecondGood)/static_cast<double>(CountMCTruthMuon) <<"%" << std::endl;
1192 
1193  std::cout << "CountBadLeadingMuonTrackAndOnlyOneMuonTrack: " << "\t" << CountBadLeadingMuonTrackAndOnlyOneMuonTrack << " = " << 100*static_cast<double>(CountBadLeadingMuonTrackAndOnlyOneMuonTrack)/static_cast<double>(CountMCTruthMuon) <<"%" << std::endl;
1194 
1195 
1196  std::cout << "CountBadLeadingMuonTrackAndLeadingPlusSecondBad: " << "\t" << CountBadLeadingMuonTrackAndLeadingPlusSecondBad << " = " << 100*static_cast<double>(CountBadLeadingMuonTrackAndLeadingPlusSecondBad)/static_cast<double>(CountMCTruthMuon) <<"%" << std::endl;
1197 
1198  std::cout << "CountNoRecoTracks: " << "\t" << CountNoRecoTracks << " = " << 100*static_cast<double>(CountNoRecoTracks)/static_cast<double>(CountMCTruthMuon) <<"%" << std::endl;
1199 
1200  std::cout << "CountNoMuonTracks: " << "\t" << CountNoMuonTracks << " = " << 100*static_cast<double>(CountNoMuonTracks)/static_cast<double>(CountMCTruthMuon) <<"%" << std::endl;
1201 
1202  std::cout << std::endl;
1203 
1204  std::cout << "CountBadLeadingMuonTrackAndOnlyOneMuonTrackCompleteness: " << CountBadLeadingMuonTrackAndOnlyOneMuonTrackCompleteness << std::endl;
1205  std::cout << "CountBadLeadingMuonTrackAndOnlyOneMuonTrackPurity: " << CountBadLeadingMuonTrackAndOnlyOneMuonTrackPurity << std::endl;
1206  std::cout << "CountBadLeadingMuonTrackAndOnlyOneMuonTrackTrackTooShort: " << CountBadLeadingMuonTrackAndOnlyOneMuonTrackTrackTooShort << std::endl;
1207  std::cout << "CountBadLeadingMuonTrackAndOnlyOneMuonTrackTrackTooLong: " << CountBadLeadingMuonTrackAndOnlyOneMuonTrackTrackTooLong << std::endl;
1208 
1209  std::cout << std::endl;
1210 
1211  std::cout << "CountBadLeadingMuonTrackAndLeadingPlusSecondBadCompleteness: " << "\t" << CountBadLeadingMuonTrackAndLeadingPlusSecondBadCompleteness << " = " << 100*static_cast<double>(CountBadLeadingMuonTrackAndLeadingPlusSecondBadCompleteness)/static_cast<double>(CountMCTruthMuon) <<"%" << std::endl;
1212 
1213  std::cout << "CountBadLeadingMuonTrackAndLeadingPlusSecondBadPurity: " << "\t" << CountBadLeadingMuonTrackAndLeadingPlusSecondBadPurity << " = " << 100*static_cast<double>(CountBadLeadingMuonTrackAndLeadingPlusSecondBadPurity)/static_cast<double>(CountMCTruthMuon) <<"%" << std::endl;
1214 
1215  std::cout << "CountBadLeadingMuonTrackAndLeadingPlusSecondBadTrackTooShort: " << "\t" << CountBadLeadingMuonTrackAndLeadingPlusSecondBadTrackTooShort << " = " << 100*static_cast<double>(CountBadLeadingMuonTrackAndLeadingPlusSecondBadTrackTooShort)/static_cast<double>(CountMCTruthMuon) <<"%" << std::endl;
1216 
1217  std::cout << "CountBadLeadingMuonTrackAndLeadingPlusSecondBadTrackTooLong: " << "\t" << CountBadLeadingMuonTrackAndLeadingPlusSecondBadTrackTooLong << " = " << 100*static_cast<double>(CountBadLeadingMuonTrackAndLeadingPlusSecondBadTrackTooLong)/static_cast<double>(CountMCTruthMuon) <<"%" << std::endl;
1218 
1219  std::cout << std::endl;
1220 
1221  std::cout << "GoodEvents1MuonTrack: " << GoodEvents1MuonTrack << std::endl;
1222  std::cout << "GoodEvents2MuonTrack: " << GoodEvents2MuonTrack << std::endl;
1223  std::cout << "GoodEvents3MuonTrack: " << GoodEvents3MuonTrack << std::endl;
1224  std::cout << "GoodEvents4OrMoreMuonTrack: " << GoodEvents4OrMoreMuonTrack << std::endl;
1225 
1226  std::cout << "BadEvents0MuonTrack: " << BadEvents0MuonTrack << std::endl;
1227  std::cout << "BadEvents1MuonTrack: " << BadEvents1MuonTrack << std::endl;
1228  std::cout << "BadEvents2MuonTrack: " << BadEvents2MuonTrack << std::endl;
1229  std::cout << "BadEvents3MuonTrack: " << BadEvents3MuonTrack << std::endl;
1230  std::cout << "BadEvents4OrMoreMuonTrack: " << BadEvents4OrMoreMuonTrack << std::endl;
1231 
1232  art::ServiceHandle<art::TFileService> tfs;
1233 
1234  for(int i=0; i < (NCriteriaBins+1)/2; ++i)
1235  {
1236  for(int j=0; j < (NRecoTracksBins+1)/2; ++j)
1237  {
1238  if(i==0)
1239  {
1240  h_Criteria_NRecoTrack_den->SetBinContent(1+2*i,1+2*j,CountNoRecoTracks);
1241  h_Criteria_NMuonTrack_den->SetBinContent(1+2*i,1+2*j,CountNoRecoTracks);
1242  }
1243  else if(i==1)
1244  {
1245  h_Criteria_NRecoTrack_den->SetBinContent(1+2*i,1+2*j,CountNoMuonTracks);
1246  h_Criteria_NMuonTrack_den->SetBinContent(1+2*i,1+2*j,CountNoMuonTracks);
1247  }
1248  else if(i==2)
1249  {
1250  h_Criteria_NRecoTrack_den->SetBinContent(1+2*i,1+2*j,CountCompleteness);
1251  h_Criteria_NMuonTrack_den->SetBinContent(1+2*i,1+2*j,CountCompleteness);
1252  }
1253  else if(i==3)
1254  {
1255  h_Criteria_NRecoTrack_den->SetBinContent(1+2*i,1+2*j,CountPurity);
1256  h_Criteria_NMuonTrack_den->SetBinContent(1+2*i,1+2*j,CountPurity);
1257  }
1258  else if(i==4)
1259  {
1260  h_Criteria_NRecoTrack_den->SetBinContent(1+2*i,1+2*j,CountTrackLengthTooShort);
1261  h_Criteria_NMuonTrack_den->SetBinContent(1+2*i,1+2*j,CountTrackLengthTooShort);
1262  }
1263  else if(i==5)
1264  {
1265  h_Criteria_NRecoTrack_den->SetBinContent(1+2*i,1+2*j,CountTrackLengthTooLong);
1266  h_Criteria_NMuonTrack_den->SetBinContent(1+2*i,1+2*j,CountTrackLengthTooLong);
1267  }
1268  else if(i==6)
1269  {
1270  h_Criteria_NRecoTrack_den->SetBinContent(1+2*i,1+2*j,CountRecoMuon);
1271  h_Criteria_NMuonTrack_den->SetBinContent(1+2*i,1+2*j,CountRecoMuon);
1272  }
1273  }
1274  }
1275 
1276 
1277  h_Efficiency_ThetaXZ->Divide(h_ThetaXZ_num,h_ThetaXZ_den,1,1, "");
1278  h_Efficiency_ThetaYZ->Divide(h_ThetaYZ_num,h_ThetaYZ_den,1,1, "");
1279  h_Efficiency_SinThetaYZ->Divide(h_SinThetaYZ_num,h_SinThetaYZ_den,1,1, "");
1280 
1281  h_Efficiency_ThetaXZ_ThetaYZ->Divide(h_ThetaXZ_ThetaYZ_num, h_ThetaXZ_ThetaYZ_den, 1, 1, "");
1282  h_Efficiency_ThetaXZ_SinThetaYZ->Divide(h_ThetaXZ_SinThetaYZ_num, h_ThetaXZ_SinThetaYZ_den, 1, 1, "");
1283 
1284  h_Criteria_NRecoTrack->Divide(h_Criteria_NRecoTrack_num, h_Criteria_NRecoTrack_den, 1, 1, "");
1285  h_Criteria_NMuonTrack->Divide(h_Criteria_NMuonTrack_num, h_Criteria_NMuonTrack_den, 1, 1, "");
1286 
1287  h_ThetaXZ_ThetaYZ_LeadingPlusSecondOk_num->Add(h_ThetaXZ_ThetaYZ_num,1);
1288  h_ThetaXZ_ThetaYZ_LeadingPlusSecondOk_num->Add(h_ThetaXZ_ThetaYZ_LeadingPlusSecondOk,1);
1289  h_Efficiency_ThetaXZ_ThetaYZ_LeadingPlusSecond->Divide(h_ThetaXZ_ThetaYZ_LeadingPlusSecondOk_num,h_ThetaXZ_ThetaYZ_den, 1, 1, "");
1290 
1291  h_ThetaXZ_SinThetaYZ_LeadingPlusSecondOk_num->Add(h_ThetaXZ_SinThetaYZ_num,1);
1292  h_ThetaXZ_SinThetaYZ_LeadingPlusSecondOk_num->Add(h_ThetaXZ_SinThetaYZ_LeadingPlusSecondOk,1);
1293  h_Efficiency_ThetaXZ_SinThetaYZ_LeadingPlusSecond->Divide(h_ThetaXZ_SinThetaYZ_LeadingPlusSecondOk_num,h_ThetaXZ_SinThetaYZ_den, 1, 1, "");
1294 
1295  h_Efficiency_ThetaXZ_ThetaYZ_DifferenceLeadingAndLeadingPlusSecond->Add(h_Efficiency_ThetaXZ_ThetaYZ_LeadingPlusSecond,1);
1296  h_Efficiency_ThetaXZ_ThetaYZ_DifferenceLeadingAndLeadingPlusSecond->Add(h_Efficiency_ThetaXZ_ThetaYZ,-1);
1297 }
1298 //========================================================================
1299 DEFINE_ART_MODULE(MuonTrackingEff)
1300 
1301 }
1302 
1303 #endif // MuonTrackingEff_Module
x
Float_t x
Definition: compare.C:6
LArG4Parameters.h
Store parameters for running LArG4.
cheat::ParticleInventoryService::TrackIdToParticle_P
const simb::MCParticle * TrackIdToParticle_P(int const &id)
Definition: ParticleInventoryService_service.cc:134
simb::MCParticle::NumberTrajectoryPoints
unsigned int NumberTrajectoryPoints() const
Definition: MCParticle.h:222
simb::MCParticle::Position
const TLorentzVector & Position(const int i=0) const
Definition: MCParticle.h:223
simb::MCParticle::PdgCode
int PdgCode() const
Definition: MCParticle.h:216
detinfo::DetectorClocks::TPCClock
virtual const ::detinfo::ElecClock & TPCClock() const =0
Lends a constant TPC clock with time set to trigger time.
DUNE::MuonTrackingEff::FuncDistanceAndAngleBetweenTruthAndRecoTrack
void FuncDistanceAndAngleBetweenTruthAndRecoTrack(const simb::MCParticle *&MCparticle, art::Ptr< recob::Track > Track, double &TempDistanceBetweenTruthAndRecoTrack, double &TempAngleBeetweenTruthAndRecoTrack)
Definition: MuonTrackingEff_module.cc:1094
recob::Track::VertexDirection
TVector3 VertexDirection() const
Covariance matrices are either set or not.
Definition: Track.h:247
iterator
intermediate_table::iterator iterator
Definition: intermediate_table.cc:24
cheat::BackTrackerService::HitToTrackIDEs
const std::vector< sim::TrackIDE > HitToTrackIDEs(recob::Hit const &hit)
Definition: BackTrackerService_service.cc:145
mf::LogInfo
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
Definition: MessageLogger.h:212
TrackIDs
std::vector< TrackID > TrackIDs
Definition: ClusteringValidation_module.cc:75
simb::MCParticle::Mother
int Mother() const
Definition: MCParticle.h:217
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:279
DUNE::MuonTrackingEff::h_MuonTrackStitching_FailedCriteria_CriteriaTwoTracks_Angle
TH2D * h_MuonTrackStitching_FailedCriteria_CriteriaTwoTracks_Angle
Definition: MuonTrackingEff_module.cc:242
geo::TPCGeo
Geometry information for a single TPC.
Definition: TPCGeo.h:37
DUNE::MuonTrackingEff::h_MuonTrackStitching_MatchedCriteria_TrackResLeading_TrackResSecond
TH2D * h_MuonTrackStitching_MatchedCriteria_TrackResLeading_TrackResSecond
Definition: MuonTrackingEff_module.cc:246
DUNE::MuonTrackingEff::h_Efficiency_ThetaXZ_SinThetaYZ_LeadingPlusSecond
TH2D * h_Efficiency_ThetaXZ_SinThetaYZ_LeadingPlusSecond
Definition: MuonTrackingEff_module.cc:202
DUNE::MuonTrackingEff::h_MuonTrackStitching_FailedCriteria_Distance_Angle
TH2D * h_MuonTrackStitching_FailedCriteria_Distance_Angle
Definition: MuonTrackingEff_module.cc:239
art::DataViewImpl::get
bool get(SelectorBase const &, Handle< PROD > &result) const
Definition: DataViewImpl.h:307
std
STL namespace.
DUNE::MuonTrackingEff::h_MuonTrackStitching_MatchedCriteria_TrackRes_Completeness
TH2D * h_MuonTrackStitching_MatchedCriteria_TrackRes_Completeness
Definition: MuonTrackingEff_module.cc:245
DUNE::MuonTrackingEff::reconfigure
void reconfigure(fhicl::ParameterSet const &pset)
Definition: MuonTrackingEff_module.cc:313
art::Event::isRealData
bool isRealData() const
Definition: Event.h:83
DUNE::MuonTrackingEff::geom
art::ServiceHandle< geo::Geometry > geom
Definition: MuonTrackingEff_module.cc:276
MCParticle.h
Particle class.
art::Run
Definition: Run.h:30
simb::MCParticle::TrackId
int TrackId() const
Definition: MCParticle.h:214
DUNE::MuonTrackingEff::h_MuonTrackStitching_MatchedCriteria_CriteriaTwoTracks_Angle
TH2D * h_MuonTrackStitching_MatchedCriteria_CriteriaTwoTracks_Angle
Definition: MuonTrackingEff_module.cc:248
DUNE::MuonTrackingEff::h_MuonTrackStitching_FailedCriteria_CompletenessSecondMuon_Angle
TH2D * h_MuonTrackStitching_FailedCriteria_CompletenessSecondMuon_Angle
Definition: MuonTrackingEff_module.cc:241
DUNE::MuonTrackingEff::h_MuonTrackStitching_FailedCriteria_CompletenessLeading_CompletenessSecond
TH2D * h_MuonTrackStitching_FailedCriteria_CompletenessLeading_CompletenessSecond
Definition: MuonTrackingEff_module.cc:238
DUNE::MuonTrackingEff::insideFV
bool insideFV(double vertex[4])
Definition: MuonTrackingEff_module.cc:1146
Event.h
geo::GeometryCore::FindTPCAtPosition
geo::TPCID FindTPCAtPosition(double const worldLoc[3]) const
Returns the ID of the TPC at specified location.
Definition: GeometryCore.h:2417
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:236
DUNE::MuonTrackingEff::h_MuonTrackStitching_TrackResLeading_TrackResSecond
TH2D * h_MuonTrackStitching_TrackResLeading_TrackResSecond
Definition: MuonTrackingEff_module.cc:229
recob::Track::Length
double Length(size_t p=0) const
Access to various track properties.
Definition: Track.h:174
DUNE::MuonTrackingEff::h_MuonTrackStitching_FailedCriteriaAndLeadingPlusSecondBad_TrackResLeading_TrackResSecond
TH2D * h_MuonTrackStitching_FailedCriteriaAndLeadingPlusSecondBad_TrackResLeading_TrackResSecond
Definition: MuonTrackingEff_module.cc:257
DEFINE_ART_MODULE
#define DEFINE_ART_MODULE(klass)
Definition: ModuleMacros.h:42
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:199
fhicl::ParameterSet::get
T get(std::string const &key) const
Definition: ParameterSet.h:231
DUNE::MuonTrackingEff::CountBadLeadingMuonTrackAndLeadingPlusSecondBadTrackTooShort
int CountBadLeadingMuonTrackAndLeadingPlusSecondBadTrackTooShort
Definition: MuonTrackingEff_module.cc:132
DUNE::MuonTrackingEff::h_MuonTrackStitching_FailedCriteriaAndLeadingPlusSecondBad_TrackRes_Completeness
TH2D * h_MuonTrackStitching_FailedCriteriaAndLeadingPlusSecondBad_TrackRes_Completeness
Definition: MuonTrackingEff_module.cc:256
DUNE::MuonTrackingEff::processEff
void processEff(const art::Event &evt, bool &isFiducial)
Definition: MuonTrackingEff_module.cc:613
sim::ParticleList::begin
iterator begin()
Definition: ParticleList.h:305
TrackID
TrackID
Definition: ClusteringValidation_module.cc:73
detinfo::DetectorProperties::NumberTimeSamples
virtual unsigned int NumberTimeSamples() const =0
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:1064
geo::TPCID
The data type to uniquely identify a TPC.
Definition: geo_types.h:195
Track.h
Provides recob::Track data product.
geo_types.h
Definition of data types for geometry description.
hit
Detector simulation of raw signals on wires.
Definition: DumpHits_module.cc:25
detinfo::DetectorClocks
Conversion of times between different formats and references.
Definition: DetectorClocks.h:298
DUNE::MuonTrackingEff::h_MuonTrackStitching_FailedCriteria_TrackResSecondMuon_Angle
TH2D * h_MuonTrackStitching_FailedCriteria_TrackResSecondMuon_Angle
Definition: MuonTrackingEff_module.cc:240
simb::MCParticle::Vx
double Vx(const int i=0) const
Definition: MCParticle.h:225
recob::Track::ID
int ID() const
Definition: Track.h:205
art::TFileDirectory::make
T * make(ARGS...args) const
Definition: TFileDirectory.h:82
DUNE::MuonTrackingEff::h_MuonTrackStitching_FailedCriteriaButLeadingPlusSecondGood_TrackResLeading_TrackResSecond
TH2D * h_MuonTrackStitching_FailedCriteriaButLeadingPlusSecondGood_TrackResLeading_TrackResSecond
Definition: MuonTrackingEff_module.cc:252
detinfo::ElecClock::TickPeriod
double TickPeriod() const
A single tick period in microseconds.
Definition: ElecClock.h:342
DUNE::MuonTrackingEff::CountBadLeadingMuonTrackAndLeadingPlusSecondBadTrackTooLong
int CountBadLeadingMuonTrackAndLeadingPlusSecondBadTrackTooLong
Definition: MuonTrackingEff_module.cc:133
art::DataViewImpl::getByLabel
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
detinfo::DetectorProperties::DriftVelocity
virtual double DriftVelocity(double efield=0., double temperature=0.) const =0
DUNE::MuonTrackingEff::h_MuonTrackStitching_FailedCriteriaAndLeadingPlusSecondBad_Distance_Angle
TH2D * h_MuonTrackStitching_FailedCriteriaAndLeadingPlusSecondBad_Distance_Angle
Definition: MuonTrackingEff_module.cc:258
recob::Track::Vertex
TVector3 Vertex() const
Covariance matrices are either set or not.
Definition: Track.h:245
simb::MCParticle::Momentum
const TLorentzVector & Momentum(const int i=0) const
Definition: MCParticle.h:224
DUNE::MuonTrackingEff::beginRun
void beginRun(const art::Run &run)
Definition: MuonTrackingEff_module.cc:602
DUNE::MuonTrackingEff::h_MuonTrackStitching_CompletenessSecondMuon_Angle
TH2D * h_MuonTrackStitching_CompletenessSecondMuon_Angle
Definition: MuonTrackingEff_module.cc:232
recob::Track::EndDirection
TVector3 EndDirection() const
Covariance matrices are either set or not.
Definition: Track.h:248
simb::MCParticle::Vz
double Vz(const int i=0) const
Definition: MCParticle.h:227
LOG_DEBUG
#define LOG_DEBUG(id)
Definition: MessageLogger.h:263
DUNE::MuonTrackingEff::h_MuonTrackStitching_FailedCriteriaButLeadingPlusSecondGood_TrackRes_Completeness
TH2D * h_MuonTrackStitching_FailedCriteriaButLeadingPlusSecondGood_TrackRes_Completeness
Definition: MuonTrackingEff_module.cc:251
DUNE::MuonTrackingEff::truthLength
double truthLength(const simb::MCParticle *MCparticle)
Definition: MuonTrackingEff_module.cc:1103
DUNE::MuonTrackingEff::h_MuonTrackStitching_FailedCriteria_TrackResLeading_TrackResSecond
TH2D * h_MuonTrackStitching_FailedCriteria_TrackResLeading_TrackResSecond
Definition: MuonTrackingEff_module.cc:237
DUNE::MuonTrackingEff::truthMatcher
void truthMatcher(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:1008
art::fill_ptr_vector
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
Definition: Ptr.h:464
e
Float_t e
Definition: plot.C:34
geo
Namespace collecting geometry-related classes utilities.
Definition: TrackContainmentAlg.hh:20
track
Float_t track
Definition: plot.C:34
recob::Track::End
TVector3 End() const
Covariance matrices are either set or not.
Definition: Track.h:246
Track
Definition: Track.hh:43
DUNE::MuonTrackingEff::CountBadLeadingMuonTrackAndLeadingPlusSecondBadCompleteness
int CountBadLeadingMuonTrackAndLeadingPlusSecondBadCompleteness
Definition: MuonTrackingEff_module.cc:130
geo::TPCGeo::LocalToWorld
void LocalToWorld(const double *tpc, double *world) const
Transform point from local TPC frame to world frame.
Definition: TPCGeo.h:490
simb::MCParticle::Vy
double Vy(const int i=0) const
Definition: MCParticle.h:226
DUNE::MuonTrackingEff::h_MuonTrackStitching_MatchedCriteria_Distance_Angle
TH2D * h_MuonTrackStitching_MatchedCriteria_Distance_Angle
Definition: MuonTrackingEff_module.cc:247
event
Event finding and building.
Definition: AggregateEvent_module.cc:55
DUNE::MuonTrackingEff::analyze
void analyze(const art::Event &evt)
Definition: MuonTrackingEff_module.cc:606
DUNE::MuonTrackingEff::h_MuonTrackStitching_FailedCriteriaButLeadingPlusSecondGood_Distance_Angle
TH2D * h_MuonTrackStitching_FailedCriteriaButLeadingPlusSecondGood_Distance_Angle
Definition: MuonTrackingEff_module.cc:253
vertex
vertex reconstruction
Definition: AggregateVertex_module.cc:43