30 #include "art_root_io/TFileService.h" 31 #include "canvas/Persistency/Common/FindManyP.h" 41 #define MAX_TRACKS 1000 54 void endJob()
override;
58 void processEff(
const art::Event& evt,
bool& isFiducial);
66 double& TotalRecoEnergy);
70 double& TempDistanceBetweenTracks,
71 double& TempAngleBetweenTracks,
72 double& TempCriteriaTwoTracks);
74 void FuncDistanceAndAngleBetweenTruthAndRecoTrack(
const simb::MCParticle*& MCparticle,
76 double& TempDistanceBetweenTruthAndRecoTrack,
77 double& TempAngleBeetweenTruthAndRecoTrack);
81 bool insideFV(
double vertex[4]);
83 void doEfficiencies();
93 double MCTruthMuonVertex[4];
101 double MCTruthMuonThetaXZ = 0;
102 double MCTruthMuonThetaYZ = 0;
105 int EventCounter = 0;
107 int CountMCTruthMuon = 0;
108 int CountRecoMuon = 0;
110 int CountGoodLeadingMuonTrack = 0;
111 int CountNoRecoTracks = 0;
112 int CountNoMuonTracks = 0;
113 int CountBadLeadingMuonTrack = 0;
114 int CountCompleteness = 0;
116 int CountTrackLengthTooShort = 0;
117 int CountTrackLengthTooLong = 0;
119 int CountBadLeadingMuonTrackAndOnlyOneMuonTrack = 0;
120 int CountBadLeadingMuonTrackButLeadingPlusSecondGood = 0;
121 int CountBadLeadingMuonTrackAndLeadingPlusSecondBad = 0;
123 int CountBadLeadingMuonTrackAndLeadingPlusSecondBadCompleteness = 0;
124 int CountBadLeadingMuonTrackAndLeadingPlusSecondBadPurity = 0;
125 int CountBadLeadingMuonTrackAndLeadingPlusSecondBadTrackTooShort = 0;
126 int CountBadLeadingMuonTrackAndLeadingPlusSecondBadTrackTooLong = 0;
128 int CountBadLeadingMuonTrackAndOnlyOneMuonTrackCompleteness = 0;
129 int CountBadLeadingMuonTrackAndOnlyOneMuonTrackPurity = 0;
130 int CountBadLeadingMuonTrackAndOnlyOneMuonTrackTrackTooShort = 0;
131 int CountBadLeadingMuonTrackAndOnlyOneMuonTrackTrackTooLong = 0;
136 int GoodEvents1MuonTrack = 0;
137 int GoodEvents2MuonTrack = 0;
138 int GoodEvents3MuonTrack = 0;
139 int GoodEvents4OrMoreMuonTrack = 0;
141 int BadEvents0MuonTrack = 0;
142 int BadEvents1MuonTrack = 0;
143 int BadEvents2MuonTrack = 0;
144 int BadEvents3MuonTrack = 0;
145 int BadEvents4OrMoreMuonTrack = 0;
267 int NThetaXZBins = 36;
268 int ThetaXZBinMin = 0;
269 int ThetaXZBinMax = 360;
271 int NThetaYZBins = 18;
272 int ThetaYZBinMin = -90;
273 int ThetaYZBinMax = 90;
275 int NSinThetaYZBins = 18;
276 int SinThetaYZBinMin = -1;
277 int SinThetaYZBinMax = 1;
279 int NCriteriaBins = 13;
280 double CriteriaBinMin = -0.25;
281 double CriteriaBinMax = 6.25;
283 int NRecoTracksBins = 19;
284 double RecoTracksBinMin = -0.25;
285 double RecoTracksBinMax = 9.25;
303 std::cout <<
"job begin..." <<
std::endl;
305 auto const*
geo = lar::providerFrom<geo::Geometry>();
314 for (
size_t i = 0; i <
geo->NTPC(); ++i) {
315 double local[3] = {0., 0., 0.};
316 double world[3] = {0., 0., 0.};
319 if (minx > world[0] -
geo->DetHalfWidth(i)) minx = world[0] -
geo->DetHalfWidth(i);
320 if (maxx < world[0] +
geo->DetHalfWidth(i)) maxx = world[0] +
geo->DetHalfWidth(i);
321 if (miny > world[1] -
geo->DetHalfHeight(i)) miny = world[1] -
geo->DetHalfHeight(i);
322 if (maxy < world[1] +
geo->DetHalfHeight(i)) maxy = world[1] +
geo->DetHalfHeight(i);
323 if (minz > world[2] -
geo->DetLength(i) / 2.) minz = world[2] -
geo->DetLength(i) / 2.;
324 if (maxz < world[2] +
geo->DetLength(i) / 2.) maxz = world[2] +
geo->DetLength(i) / 2.;
334 std::cout <<
"Fiducial volume:" 343 gStyle->SetTitleOffset(1.3,
"Y");
347 tfs->make<TH1D>(
"h_Purity",
"All events: Purity vs. # events; Purity; # events", 60, 0, 1.2);
350 "h_Completeness",
"All events: Completeness vs # events; Completeness; # events", 60, 0, 1.2);
354 tfs->make<TH1D>(
"h_TrackRes",
355 "All events: L_{reco}/L_{truth} vs. # events; L_{reco}/L_{truth}; # events;",
362 "All events: Total reco energy (sum of all hits in all " 363 "tracks) vs. # events; E_{reco., tot.} [MeV]; # events",
370 "All events: truth muon length vs. # events; truth muon length [cm]; # events",
377 "All events: Vertex residuals vs. # events; #Delta vertex_{truth-teco} [cm]; # events",
384 "All events: Angular residuals vs. # events; #Delta#theta_{truth-reco} [#circ]; # events",
391 tfs->make<TH1D>(
"h_Efficiency_ThetaXZ",
392 "Muon reco efficiency vs. #theta_{XZ}; #theta_{XZ} [#circ]; Efficiency",
397 tfs->make<TH1D>(
"h_ThetaXZ_den",
398 "# generated muons vs. #theta_{XZ}; #theta_{XZ} [#circ]; # generated muons",
404 "# reconstructed muons vs. #theta_{XZ}; #theta_{XZ} [#circ]; # reconstructed muons",
411 "h_Efficiency_ThetaYZ",
412 "Muon reco efficiency vs. #theta_{YZ}; #theta_{YZ} [#circ]; Muon reco efficiency",
418 tfs->make<TH1D>(
"h_ThetaYZ_den",
419 "# generated muons vs. #theta_{YZ}; #theta_{YZ} [#circ]; # generated muons",
425 "# reconstructed muons vs. #theta_{YZ}; #theta_{YZ} [#circ]; # reconstructed muons",
432 "h_Efficiency_SinThetaYZ",
433 "Muon reco efficiency vs. sin(#theta_{YZ}); sin(#theta_{YZ}); Muon reco efficiency",
439 tfs->make<TH1D>(
"h_SinThetaYZ_den",
440 "# generated muons vs. sin(#theta_{YZ}); sin(#theta_{YZ}); # generated muons",
446 "# reconstructed muons vs. sin(#theta_{YZ}); sin(#theta_{YZ}); # reconstructed muons",
455 h_TruthLength->Sumw2();
456 h_VertexRes->Sumw2();
457 h_DirectionRes->Sumw2();
459 h_Efficiency_SinThetaYZ->Sumw2();
461 h_SinThetaYZ_num->Sumw2();
463 h_Efficiency_ThetaXZ->Sumw2();
464 h_ThetaXZ_den->Sumw2();
465 h_ThetaXZ_num->Sumw2();
467 h_Efficiency_ThetaYZ->Sumw2();
469 h_ThetaYZ_num->Sumw2();
474 "h_Efficiency_ThetaXZ_ThetaYZ",
475 "Muon reco efficiency: #theta_{XZ} vs. #theta_{YZ}; #theta_{XZ} [#circ]; #theta_{YZ} [#circ]",
485 "h_ThetaXZ_ThetaYZ_den",
486 "# generated muons: #theta_{XZ} vs. #theta_{YZ}; #theta_{XZ} [#circ]; #theta_{YZ} [#circ]",
496 "# reconstructed muons: #theta_{XZ} vs. #theta_{YZ}; " 497 "#theta_{XZ} [#circ]; #theta_{YZ} [#circ]",
507 tfs->make<TH2D>(
"h_FailedReconstruction_ThetaXZ_ThetaYZ",
508 "# failed reconstructions: #theta_{XZ} vs. #theta_{YZ}; #theta_{XZ} [#circ]; " 509 "#theta_{YZ} [#circ]",
520 tfs->make<TH2D>(
"h_Efficiency_ThetaXZ_SinThetaYZ",
521 "Muon reco efficiency: #theta_{XZ} vs. sin(#theta_{YZ}); #theta_{XZ} " 522 "[#circ]; sin(#theta_{YZ})",
532 "h_ThetaXZ_SinThetaYZ_den",
533 "# generated muons: #theta_{XZ} vs. sin(#theta_{YZ}); #theta_{XZ} [#circ]; sin(#theta_{YZ})",
543 tfs->make<TH2D>(
"h_ThetaXZ_SinThetaYZ_num",
544 "# reconstructed muons: #theta_{XZ} vs. sin(#theta_{YZ}); #theta_{XZ} " 545 "[#circ]; sin(#theta_{YZ}); #theta_{XZ} [#circ]; sin(#theta_{YZ})",
555 tfs->make<TH2D>(
"h_FailedReconstruction_ThetaXZ_SinThetaYZ",
556 "# failed reconstructions: #theta_{XZ} vs. sin(#theta_{YZ}); #theta_{XZ} " 557 "[#circ]; sin(#theta_{YZ})",
568 tfs->make<TH2D>(
"h_Efficiency_ThetaXZ_ThetaYZ_LeadingPlusSecond",
569 "Muon reco efficiency after stitching: #theta_{XZ} vs. #theta_{YZ}; " 570 "#theta_{XZ} [#circ]; #theta_{YZ} [#circ]",
580 tfs->make<TH2D>(
"h_ThetaXZ_ThetaYZ_LeadingPlusSecondOk",
581 "# reconstructed muons after stitching (failed before stitching): " 582 "#theta_{XZ} vs #theta_{YZ}; #theta_{XZ} [#circ]; #theta_{YZ} [#circ]",
592 tfs->make<TH2D>(
"h_ThetaXZ_ThetaYZ_LeadingPlusSecondOk_num",
593 "# reconstructed muons after stitching: #theta_{XZ} vs. #theta_{YZ}; " 594 "#theta_{XZ} [#circ]; #theta_{YZ} [#circ]",
605 tfs->make<TH2D>(
"h_Efficiency_ThetaXZ_SinThetaYZ_LeadingPlusSecond",
606 "Muon reco efficiency after stitching: #theta_{XZ} vs. sin(#theta_{YZ}); " 607 "#theta_{XZ} [#circ]; sin(#theta_{YZ})",
617 tfs->make<TH2D>(
"h_ThetaXZ_SinThetaYZ_LeadingPlusSecondOk",
618 "# reconstructed muons after stitching (failed before stitching): " 619 "#theta_{XZ} vs. sin(#theta_{YZ}); #theta_{XZ} [#circ]; sin(#theta_{YZ})",
629 tfs->make<TH2D>(
"h_ThetaXZ_SinThetaYZ_LeadingPlusSecondOk_num",
630 "# reconstructed muons after stitching: #theta_{XZ} vs. sin(#theta_{YZ}); " 631 "#theta_{XZ} [#circ]; sin(#theta_{YZ})",
642 tfs->make<TH2D>(
"h_Efficiency_ThetaXZ_ThetaYZ_DifferenceLeadingAndLeadingPlusSecond",
643 "Muon reco efficiency: difference before and after stitching: #theta_{XZ} " 644 "vs. #theta_{YZ}; #theta_{XZ} [#circ]; #theta_{YZ} [#circ]",
655 tfs->make<TH2D>(
"h_NoRecoTrackAtAll_ThetaXZ_SinThetaYZ",
656 "# events with no reco track at all: #theta_{XZ} vs. sin(#theta_{YZ}); " 657 "#theta_{XZ} [#circ]; sin(#theta_{YZ})",
667 tfs->make<TH2D>(
"h_NoMuonTrack_ThetaXZ_SinThetaYZ",
668 "# events with no muon track: #theta_{XZ} vs. sin(#theta_{YZ}); #theta_{XZ} " 669 "[#circ]; sin(#theta_{YZ})",
679 tfs->make<TH2D>(
"h_TrackTooShort_ThetaXZ_SinThetaYZ",
680 "# events with L_{reco}/L_{truth} < 75%: #theta_{XZ} vs. sin(#theta_{YZ}); " 681 "#theta_{XZ} [#circ]; sin(#theta_{YZ})",
691 tfs->make<TH2D>(
"h_TrackTooLong_ThetaXZ_SinThetaYZ",
692 "#events with L_{reco}/L_{truth} > 125%: #theta_{XZ} vs. sin(#theta_{YZ}); " 693 "#theta_{XZ} [#circ]; sin(#theta_{YZ})",
703 tfs->make<TH2D>(
"h_Completeness_ThetaXZ_SinThetaYZ",
704 "# events with Completeness < 50%: #theta_{XZ} vs. sin(#theta_{YZ}); " 705 "#theta_{XZ} [#circ]; sin(#theta_{YZ})",
715 tfs->make<TH2D>(
"h_Purity_ThetaXZ_SinThetaYZ",
716 "# events with Purity < 50%: #theta_{XZ} vs. sin(#theta_{YZ}); #theta_{XZ} " 717 "[#circ]; sin(#theta_{YZ})",
728 tfs->make<TH2D>(
"h_Criteria_NRecoTrack",
729 "Ratio: criteria vs. # reco tracks; Criteria; # reco tracks",
739 tfs->make<TH2D>(
"h_Criteria_NRecoTrack_num",
740 "# events: criteria vs. # reco tracks; Criteria; # reco tracks",
750 "Divider histogram; Criteria; # reco tracks",
761 tfs->make<TH2D>(
"h_Criteria_NMuonTrack",
762 "Ratio: criteria vs. # muon tracks; Criteria; # muon tracks",
772 tfs->make<TH2D>(
"h_Criteria_NMuonTrack_num",
773 "# events: criteria vs. # muon tracks; Criteria; # muon tracks",
783 "Divider histogram; Criteria; # muon tracks",
794 "h_NoMuonTrack_MaxTrackLength_PDGCode",
795 "Events with no muon track: L_{reco, max} vs. PDG Code; L_{reco} [cm]; PDG Code",
806 tfs->make<TH2D>(
"h_MuonTrackStitching_TrackRes_Completeness",
807 "All events: L_{reco}/L_{truth} (leading) vs. Completeness (leading); " 808 "L_{reco}/L_{truth} (leading); Completeness (leading)",
818 tfs->make<TH2D>(
"h_MuonTrackStitching_TrackResLeading_TrackResSecond",
819 "All events: L_{reco}/L_{truth} (leading) vs. L_{reco}/L_{truth} (second); " 820 "L_{reco}/L_{truth} (leading); L_{reco}/L_{truth} (second)",
830 tfs->make<TH2D>(
"h_MuonTrackStitching_Distance_Angle",
831 "All events: distance vs. angle b/w leading and second muon track; Distance " 832 "[cm]; Angle [#circ]",
842 tfs->make<TH2D>(
"h_MuonTrackStitching_TrackResSecondMuon_Angle",
843 "All events: L_{reco}/L_{truth} (second) vs. angle; L_{reco}/L_{truth} " 844 "(second); Angle [#circ]",
854 "h_MuonTrackStitching_CompletenessSecondMuon_Angle",
855 "All events: Completeness (second) vs. angle; Completeness (second); Angle [#circ]",
865 tfs->make<TH2D>(
"h_MuonTrackStitching_CriteriaTwoTracks_Angle",
866 "All events: CriteriaTwoTracks vs. angle; Criteria; Angle [#circ]",
877 tfs->make<TH2D>(
"h_MuonTrackStitching_FailedCriteria_TrackRes_Completeness",
878 "Bad events: L_{reco}/L_{truth} (leading) vs. Completeness (leading); " 879 "L_{reco}/L_{truth} (leading); Completeness (leading)",
889 tfs->make<TH2D>(
"h_MuonTrackStitching_FailedCriteria_Distance_Angle",
890 "Bad events: distance vs. angle b/w leading and second muon track; Distance " 891 "[cm]; Angle [#circ]",
901 tfs->make<TH2D>(
"h_MuonTrackStitching_FailedCriteria_TrackResLeading_TrackResSecond",
902 "Bad events: L_{reco}/L_{truth} (leading) vs. L_{reco}/L_{truth} (second); " 903 "L_{reco}/L_{truth} (leading); L_{reco}/L_{truth} (second)",
913 tfs->make<TH2D>(
"*h_MuonTrackStitching_FailedCriteria_CompletenessLeading_CompletenessSecond",
914 "Bad events: Completeness (leading) vs. Completeness (second); Completeness " 915 "(leading); Completeness (second)",
925 tfs->make<TH2D>(
"h_MuonTrackStitching_FailedCriteria_TrackResSecondMuon_Angle",
926 "Bad events: L_{reco}/L_{truth} (second) vs. angle; L_{reco}/L_{truth} " 927 "(second); Angle [#circ]",
937 "h_MuonTrackStitching_FailedCriteria_CompletenessSecondMuon_Angle",
938 "Bad events: Completeness (second) vs. angle; Completeness (second); Angle [#circ]",
948 tfs->make<TH2D>(
"h_MuonTrackStitching_FailedCriteria_CriteriaTwoTracks_Angle",
949 "Bad events: CriteriaTwoTracks vs. angle; Criteria; Angle [#circ]",
960 tfs->make<TH2D>(
"h_MuonTrackStitching_MatchedCriteria_TrackRes_Completeness",
961 "Good events: L_{reco}/L_{truth} (leading) vs. Completeness (leading); " 962 "L_{reco}/L_{truth} (leading); Completeness (leading)",
972 tfs->make<TH2D>(
"h_MuonTrackStitching_MatchedCriteria_Distance_Angle",
973 "Good events: distance vs. angle b/w leading and second muon track; Distance " 974 "[cm]; Angle [#circ]",
984 tfs->make<TH2D>(
"h_MuonTrackStitching_MatchedCriteria_TrackResLeading_TrackResSecond",
985 "Good events: L_{reco}/L_{truth} (leading) vs. L_{reco}/L_{truth} (second); " 986 "L_{reco}/L_{truth} (leading); L_{reco}/L_{truth} (second)",
996 tfs->make<TH2D>(
"h_MuonTrackStitching_MatchedCriteria_CriteriaTwoTracks_Angle",
997 "Good events: CriteriaTwoTracks vs. angle b/w leading and second muon track; " 998 "Criteria; Angle [#circ]",
1010 "h_MuonTrackStitching_FailedCriteriaButLeadingPlusSecondGood_TrackRes_Completeness",
1011 "Bad events but leading + second is good: L_{reco}/L_{truth} (leading) vs. Completeness " 1012 "(leading); L_{reco}/L_{truth} (leading); Completeness (leading)",
1023 tfs->make<TH2D>(
"h_MuonTrackStitching_FailedCriteriaButLeadingPlusSecondGood_Distance_Angle",
1024 "Bad events but leading + second is good: distance vs. angle b/w leading and " 1025 "second muon track; Distance [cm]; Angle [#circ]",
1036 "h_MuonTrackStitching_FailedCriteriaButLeadingPlusSecondGood_TrackResLeading_" 1038 "Bad events but leading + second is good: L_{reco}/L_{truth} (leading) vs. " 1039 "L_{reco}/L_{truth} (second); L_{reco}/L_{truth} (leading); L_{reco}/L_{truth} (second)",
1047 ->SetOption(
"colz");
1052 "h_MuonTrackStitching_FailedCriteriaAndLeadingPlusSecondBad_TrackRes_Completeness",
1053 "Bad events and leading + second is bad: L_{reco}/L_{truth} (leading) vs. Completeness " 1054 "(leading); L_{reco}/L_{truth} (leading); Completeness (leading)",
1065 tfs->make<TH2D>(
"h_MuonTrackStitching_FailedCriteriaAndLeadingPlusSecondBad_Distance_Angle",
1066 "Bad events and leading + second is bad: distance vs. angle b/w leading and " 1067 "second muon track; Distance [cm]; Angle [#circ]",
1078 "h_MuonTrackStitching_FailedCriteriaAndLeadingPlusSecondBad_TrackResLeading_TrackResSecond",
1079 "Bad events and leading + second is bad: L_{reco}/L_{truth} (leading) vs. " 1080 "L_{reco}/L_{truth} (second); L_{reco}/L_{truth} (leading); L_{reco}/L_{truth} (second)",
1088 ->SetOption(
"colz");
1109 bool isFiducial =
false;
1121 const sim::ParticleList& plist = pi_serv->
ParticleList();
1125 particle = ipar->second;
1127 if (particle->
Mother() ==
1129 const TLorentzVector& positionStart = particle->
Position(0);
1130 positionStart.GetXYZT(
1135 MCTruthMuonParticle = particle;
1143 (180.0 / 3.14159) * atan(particle->
Momentum().Px() / particle->
Momentum().Pz());
1145 else if (particle->
Momentum().Pz() < 0 && particle->
Momentum().Px() >= 0) {
1147 180.0 + (180.0 / 3.14159) * atan(particle->
Momentum().Px() / particle->
Momentum().Pz());
1151 180.0 + (180.0 / 3.14159) * atan(particle->
Momentum().Px() / particle->
Momentum().Pz());
1153 else if (particle->
Momentum().Pz() >= 0 && particle->
Momentum().Px() < 0) {
1155 360.0 + (180.0 / 3.14159) * atan(particle->
Momentum().Px() / particle->
Momentum().Pz());
1162 double MCTruthLengthMuon =
truthLength(MCTruthMuonParticle);
1169 if (!isFiducial)
return;
1172 if (MCTruthMuonParticle) {
1186 int NMuonTracks = 0;
1190 std::vector<art::Ptr<recob::Track>>
TrackList;
1192 int NRecoTracks = TrackList.size();
1193 art::FindManyP<recob::Hit> track_hits(TrackListHandle, event,
fTrackModuleLabel);
1194 if (NRecoTracks == 0) {
1195 MF_LOG_DEBUG(
"MuonTrackingEff") <<
"There are no reco tracks... bye";
1196 std::cout <<
"There are no reco tracks! MCTruthMuonThetaXZ: " <<
std::endl;
1209 MF_LOG_DEBUG(
"MuonTrackingEff") <<
"Found this many reco tracks " << NRecoTracks;
1213 double PurityLeadingMuon = 0.;
1214 double CompletenessLeadingMuon = 0.;
1215 double RecoLengthLeadingMuon = 0.;
1218 double RecoLengthSecondMuon = 0.;
1219 double CompletenessSecondMuon = 0.;
1220 double PuritySecondMuon = 0.;
1223 double TrackLengthMuonSum = 0.;
1224 double tmpTotalRecoEnergy = 0.;
1226 double MaxLengthNoRecoMuon = 0;
1227 int PDGCodeMaxLengthNoRecoMuon = 0;
1231 std::vector<art::Ptr<recob::Hit>> tmp_TrackHits = track_hits.at(0);
1232 std::vector<art::Ptr<recob::Hit>> AllHits;
1236 auto const clockData =
1240 for (
int i = 0; i < NRecoTracks; i++) {
1242 std::vector<art::Ptr<recob::Hit>> TrackHits = track_hits.at(i);
1243 double tmpPurity = 0.;
1244 double tmpCompleteness = 0.;
1248 clockData, AllHits, TrackHits, particle, tmpPurity, tmpCompleteness, tmpTotalRecoEnergy);
1251 std::cout <<
"ERROR: Truth matcher didn't find a particle!" <<
std::endl;
1257 MaxLengthNoRecoMuon = track->
Length();
1258 PDGCodeMaxLengthNoRecoMuon = particle->
PdgCode();
1262 tmpCompleteness > 0 && tmpPurity > 0) {
1265 TrackLengthMuonSum += track->
Length();
1267 if (NMuonTracks == 1) {
1268 CompletenessLeadingMuon = tmpCompleteness;
1269 PurityLeadingMuon = tmpPurity;
1270 RecoLengthLeadingMuon = track->
Length();
1271 TrackLeadingMuon = track;
1273 RecoMuonParticle = particle;
1276 if (NMuonTracks >= 2 && tmpCompleteness > CompletenessLeadingMuon) {
1278 CompletenessSecondMuon = CompletenessLeadingMuon;
1279 PuritySecondMuon = PurityLeadingMuon;
1280 RecoLengthSecondMuon = RecoLengthLeadingMuon;
1281 TrackSecondMuon = TrackLeadingMuon;
1283 CompletenessLeadingMuon = tmpCompleteness;
1284 PurityLeadingMuon = tmpPurity;
1285 RecoLengthLeadingMuon = track->
Length();
1286 TrackLeadingMuon = track;
1288 RecoMuonParticle = particle;
1291 else if (NMuonTracks >= 2 && tmpCompleteness < CompletenessLeadingMuon &&
1292 tmpCompleteness > CompletenessSecondMuon) {
1293 CompletenessSecondMuon = tmpCompleteness;
1294 PuritySecondMuon = tmpPurity;
1295 RecoLengthSecondMuon = track->
Length();
1296 TrackSecondMuon = track;
1306 if (RecoMuonParticle && MCTruthMuonParticle) {
1310 h_TrackRes->Fill(RecoLengthLeadingMuon / MCTruthLengthMuon);
1312 std::cout <<
"TrackLeadingMuon->Vertex().X(): " << TrackLeadingMuon->
Vertex().X()
1314 std::cout <<
"MCTruthMuonParticle->Vz(): " << MCTruthMuonParticle->
Vz() <<
std::endl;
1317 double DistanceBetweenTruthAndRecoTrack;
1318 double AngleBetweenTruthAndRecoTrack;
1321 DistanceBetweenTruthAndRecoTrack,
1322 AngleBetweenTruthAndRecoTrack);
1324 h_VertexRes->Fill(DistanceBetweenTruthAndRecoTrack);
1328 CompletenessLeadingMuon);
1329 if (NMuonTracks >= 2) {
1330 double DistanceBetweenTracks;
1331 double AngleBetweenTracks;
1332 double CriteriaTwoTracks;
1336 DistanceBetweenTracks,
1341 RecoLengthLeadingMuon / MCTruthLengthMuon, RecoLengthSecondMuon / MCTruthLengthMuon);
1344 RecoLengthSecondMuon / MCTruthLengthMuon, AngleBetweenTracks);
1346 AngleBetweenTracks);
1351 if (CompletenessLeadingMuon < 0.5) {
1360 if (PurityLeadingMuon < 0.5) {
1369 if (RecoLengthLeadingMuon / MCTruthLengthMuon < 0.75) {
1378 if (RecoLengthLeadingMuon / MCTruthLengthMuon > 1.25) {
1406 if (CompletenessLeadingMuon < 0.5 || PurityLeadingMuon < 0.5 ||
1407 RecoLengthLeadingMuon / MCTruthLengthMuon < 0.75 ||
1408 RecoLengthLeadingMuon / MCTruthLengthMuon > 1.25) {
1414 RecoLengthLeadingMuon / MCTruthLengthMuon, CompletenessLeadingMuon);
1421 if (NMuonTracks >= 2) {
1422 double AngleBetweenTracks;
1423 double DistanceBetweenTracks;
1424 double CriteriaTwoTracks;
1427 DistanceBetweenTracks,
1431 if (AngleBetweenTracks > 160.) {
1433 std::cout <<
"Angle b/w tracks: " << AngleBetweenTracks <<
std::endl;
1434 std::cout <<
"CriteriaTwoTracks: " << CriteriaTwoTracks <<
std::endl;
1435 std::cout <<
"CompletenessLeadingMuon: " << CompletenessLeadingMuon <<
std::endl;
1436 std::cout <<
"CompletenessSecondMuon: " << CompletenessSecondMuon <<
std::endl;
1437 std::cout <<
"PurityLeadingMuon: " << PurityLeadingMuon <<
std::endl;
1438 std::cout <<
"PuritySecondMuon: " << PuritySecondMuon <<
std::endl;
1439 std::cout <<
"TrackLeadingMuon: " << RecoLengthLeadingMuon / MCTruthLengthMuon
1441 std::cout <<
"TrackResSecondMuon: " << RecoLengthSecondMuon / MCTruthLengthMuon
1443 std::cout <<
"TrackID leading: " << TrackLeadingMuon->
ID() <<
std::endl;
1444 std::cout <<
"TrackID second: " << TrackSecondMuon->
ID() <<
std::endl;
1448 AngleBetweenTracks);
1450 RecoLengthLeadingMuon / MCTruthLengthMuon, RecoLengthSecondMuon / MCTruthLengthMuon);
1452 CompletenessLeadingMuon, CompletenessSecondMuon);
1454 RecoLengthSecondMuon / MCTruthLengthMuon, AngleBetweenTracks);
1456 CompletenessSecondMuon, AngleBetweenTracks);
1458 AngleBetweenTracks);
1459 if ((CompletenessLeadingMuon + CompletenessSecondMuon) >= 0.5 &&
1460 PurityLeadingMuon >= 0.5 && PuritySecondMuon >= 0.5 &&
1461 (RecoLengthLeadingMuon + RecoLengthSecondMuon) / MCTruthLengthMuon >= 0.75 &&
1462 (RecoLengthLeadingMuon + RecoLengthSecondMuon) / MCTruthLengthMuon <= 1.25) {
1469 RecoLengthLeadingMuon / MCTruthLengthMuon, CompletenessLeadingMuon);
1471 ->Fill(RecoLengthLeadingMuon / MCTruthLengthMuon,
1472 RecoLengthSecondMuon / MCTruthLengthMuon);
1474 DistanceBetweenTracks, AngleBetweenTracks);
1476 if ((CompletenessLeadingMuon + CompletenessSecondMuon) < 0.5 || PurityLeadingMuon < 0.5 ||
1477 PuritySecondMuon < 0.5 ||
1478 (RecoLengthLeadingMuon + RecoLengthSecondMuon) / MCTruthLengthMuon < 0.75 ||
1479 (RecoLengthLeadingMuon + RecoLengthSecondMuon) / MCTruthLengthMuon > 1.25) {
1483 RecoLengthLeadingMuon / MCTruthLengthMuon, CompletenessLeadingMuon);
1485 ->Fill(RecoLengthLeadingMuon / MCTruthLengthMuon,
1486 RecoLengthSecondMuon / MCTruthLengthMuon);
1488 DistanceBetweenTracks, AngleBetweenTracks);
1490 if ((CompletenessLeadingMuon + CompletenessSecondMuon) < 0.5) {
1493 if (PurityLeadingMuon < 0.5 || PuritySecondMuon < 0.5) {
1496 if ((RecoLengthLeadingMuon + RecoLengthSecondMuon) / MCTruthLengthMuon < 0.75) {
1499 if ((RecoLengthLeadingMuon + RecoLengthSecondMuon) / MCTruthLengthMuon > 1.25) {
1503 else if (NMuonTracks == 1) {
1505 if (CompletenessLeadingMuon < 0.5) {
1509 if (RecoLengthLeadingMuon / MCTruthLengthMuon < 0.75) {
1512 if (RecoLengthLeadingMuon / MCTruthLengthMuon > 1.25) {
1518 if (CompletenessLeadingMuon >= 0.5 && PurityLeadingMuon >= 0.5 &&
1519 RecoLengthLeadingMuon / MCTruthLengthMuon >= 0.75 &&
1520 RecoLengthLeadingMuon / MCTruthLengthMuon <= 1.25) {
1533 RecoLengthLeadingMuon / MCTruthLengthMuon, CompletenessLeadingMuon);
1540 if (NMuonTracks >= 2) {
1542 RecoLengthLeadingMuon / MCTruthLengthMuon, RecoLengthSecondMuon / MCTruthLengthMuon);
1544 double AngleBetweenTracks;
1545 double DistanceBetweenTracks;
1546 double CriteriaTwoTracks;
1549 DistanceBetweenTracks,
1553 AngleBetweenTracks);
1555 AngleBetweenTracks);
1560 if (!RecoMuonParticle && MCTruthMuonParticle) {
1572 static_cast<double>(PDGCodeMaxLengthNoRecoMuon));
1582 double& Completeness,
1583 double& TotalRecoEnergy)
1587 std::map<int, double> trkID_E;
1589 for (
size_t j = 0; j < track_hits.size(); ++j) {
1591 std::vector<sim::TrackIDE>
TrackIDs =
1600 for (
size_t k = 0;
k < TrackIDs.size();
k++) {
1601 trkID_E[TrackIDs[
k].trackID] +=
1608 double max_E = -999.0;
1609 double TotalEnergyTrack = 0.0;
1611 double PartialEnergyTrackID =
1617 if (!trkID_E.size()) {
1625 TotalEnergyTrack += ii->second;
1627 if ((ii->second) > max_E) {
1629 PartialEnergyTrackID = ii->second;
1631 TrackID = ii->first;
1634 if (TrackID < 0) E_em += ii->second;
1643 if (TrackID < 0) {
return; }
1646 Purity = PartialEnergyTrackID / TotalEnergyTrack;
1649 TotalRecoEnergy = 0;
1650 for (
size_t k = 0;
k < AllHits.size();
1655 for (
size_t l = 0;
l < TrackIDs.size(); ++
l) {
1656 if (TrackIDs[
l].trackID == TrackID)
1657 TotalRecoEnergy += TrackIDs[
l].energy;
1661 Completeness = PartialEnergyTrackID / TotalRecoEnergy;
1667 double& TempDistanceBetweenTracks,
1668 double& TempAngleBetweenTracks,
1669 double& TempCriteriaTwoTracks)
1672 TempDistanceBetweenTracks = sqrt(
pow(Track1->
End().X() - Track2->
Vertex().X(), 2) +
1675 TempAngleBetweenTracks = (180.0 / 3.14159) * Track1->
EndDirection<TVector3>().Angle(
1677 TempCriteriaTwoTracks = 1.;
1679 if (TempDistanceBetweenTracks > sqrt(
pow(Track1->
End().X() - Track2->
End().X(), 2) +
1680 pow(Track1->
End().Y() - Track2->
End().Y(), 2) +
1681 pow(Track1->
End().Z() - Track2->
End().Z(), 2))) {
1682 TempDistanceBetweenTracks = sqrt(
pow(Track1->
End().X() - Track2->
End().X(), 2) +
1683 pow(Track1->
End().Y() - Track2->
End().Y(), 2) +
1684 pow(Track1->
End().Z() - Track2->
End().Z(), 2));
1685 TempAngleBetweenTracks = 180. - (180.0 / 3.14159) * Track1->
EndDirection<TVector3>().Angle(
1687 TempCriteriaTwoTracks = 2.;
1690 if (TempDistanceBetweenTracks > sqrt(
pow(Track1->
Vertex().X() - Track2->
End().X(), 2) +
1693 TempDistanceBetweenTracks = sqrt(
pow(Track1->
Vertex().X() - Track2->
End().X(), 2) +
1696 TempAngleBetweenTracks = (180.0 / 3.14159) * Track1->
VertexDirection<TVector3>().Angle(
1698 TempCriteriaTwoTracks = 3.;
1701 if (TempDistanceBetweenTracks > sqrt(
pow(Track1->
Vertex().X() - Track2->
Vertex().X(), 2) +
1704 TempDistanceBetweenTracks = sqrt(
pow(Track1->
Vertex().X() - Track2->
Vertex().X(), 2) +
1707 TempAngleBetweenTracks = 180. - (180.0 / 3.14159) * Track1->
VertexDirection<TVector3>().Angle(
1709 TempCriteriaTwoTracks = 4.;
1717 double& TempDistanceBetweenTruthAndRecoTrack,
1718 double& TempAngleBeetweenTruthAndRecoTrack)
1720 TempDistanceBetweenTruthAndRecoTrack = sqrt(
pow(Track->
Vertex().X() - MCparticle->
Vx(), 2) +
1724 TempAngleBeetweenTruthAndRecoTrack = 0;
1734 if (!MCparticle)
return -999.0;
1736 std::vector<double> TPCLengthHits(numberTrajectoryPoints, 0);
1737 int FirstHit = 0, LastHit = 0;
1738 double TPCLength = 0.0;
1739 bool BeenInVolume =
false;
1741 for (
unsigned int MCHit = 0; MCHit < TPCLengthHits.size(); ++MCHit) {
1742 const TLorentzVector& tmpPosition = MCparticle->
Position(MCHit);
1743 double const tmpPosArray[] = {tmpPosition[0], tmpPosition[1], tmpPosition[2]};
1745 TPCLengthHits[MCHit] = sqrt(
pow((MCparticle->
Vx(MCHit - 1) - MCparticle->
Vx(MCHit)), 2) +
1746 pow((MCparticle->
Vy(MCHit - 1) - MCparticle->
Vy(MCHit)), 2) +
1747 pow((MCparticle->
Vz(MCHit - 1) - MCparticle->
Vz(MCHit)), 2));
1749 if (tpcid.isValid) {
1751 if (!BeenInVolume) {
1752 BeenInVolume =
true;
1757 for (
int Hit = FirstHit + 1;
Hit <= LastHit; ++
Hit)
1758 TPCLength += TPCLengthHits[
Hit];
1766 double x = vertex[0];
1767 double y = vertex[1];
1768 double z = vertex[2];
1782 std::cout <<
"EventCounter: " 1785 std::cout <<
"CountMCTruthMuon: " 1790 std::cout <<
"CountGoodLeadingMuonTrack (=good events): " 1796 std::cout <<
"CountBadLeadingMuonTrack+CountNoRecoTracks+CountNoMuonTracks (=bad events): " 1804 std::cout <<
"CountNoRecoTracks+CountNoMuonTracks: " 1810 std::cout <<
"CountTrackLengthTooShort: " 1816 std::cout <<
"CountCompleteness: " 1822 std::cout <<
"CountTrackLengthTooLong: " 1828 std::cout <<
"CountPurity: " 1835 std::cout <<
"GoodLeadingMuonTrack+CountBadLeadingMuonTrackButLeadingPlusSecondGood (=good " 1836 "events after stitching): " 1846 std::cout <<
"CountBadLeadingMuonTrack+CountNoRecoTracks+CountNoMuonTracks-" 1847 "CountBadLeadingMuonTrackButLeadingPlusSecondGood (=bad events after stitching) : " 1861 std::cout <<
"CountBadLeadingMuonTrackButLeadingPlusSecondGood: " 1867 std::cout <<
"CountBadLeadingMuonTrackAndOnlyOneMuonTrack: " 1873 std::cout <<
"CountBadLeadingMuonTrackAndLeadingPlusSecondBad: " 1879 std::cout <<
"CountNoRecoTracks: " 1885 std::cout <<
"CountNoMuonTracks: " 1893 std::cout <<
"CountBadLeadingMuonTrackAndOnlyOneMuonTrackCompleteness: " 1895 std::cout <<
"CountBadLeadingMuonTrackAndOnlyOneMuonTrackPurity: " 1897 std::cout <<
"CountBadLeadingMuonTrackAndOnlyOneMuonTrackTrackTooShort: " 1899 std::cout <<
"CountBadLeadingMuonTrackAndOnlyOneMuonTrackTrackTooLong: " 1904 std::cout <<
"CountBadLeadingMuonTrackAndLeadingPlusSecondBadCompleteness: " 1907 static_cast<double>(
1912 std::cout <<
"CountBadLeadingMuonTrackAndLeadingPlusSecondBadPurity: " 1918 std::cout <<
"CountBadLeadingMuonTrackAndLeadingPlusSecondBadTrackTooShort: " 1921 static_cast<double>(
1926 std::cout <<
"CountBadLeadingMuonTrackAndLeadingPlusSecondBadTrackTooLong: " 1929 static_cast<double>(
TH2D * h_MuonTrackStitching_TrackRes_Completeness
def analyze(root, level, gtrees, gbranches, doprint)
TH2D * h_NoRecoTrackAtAll_ThetaXZ_SinThetaYZ
TH2D * h_NoMuonTrack_ThetaXZ_SinThetaYZ
TH2D * h_Completeness_ThetaXZ_SinThetaYZ
void beginRun(const art::Run &run) override
int CountBadLeadingMuonTrackAndOnlyOneMuonTrackPurity
int CountBadLeadingMuonTrackButLeadingPlusSecondGood
unsigned int NumberTrajectoryPoints() const
const TLorentzVector & Position(const int i=0) const
int BadEvents4OrMoreMuonTrack
std::vector< sim::TrackIDE > HitToTrackIDEs(detinfo::DetectorClocksData const &clockData, recob::Hit const &hit) const
TH2D * h_TrackTooLong_ThetaXZ_SinThetaYZ
void FuncDistanceAndAngleBetweenTruthAndRecoTrack(const simb::MCParticle *&MCparticle, art::Ptr< recob::Track > Track, double &TempDistanceBetweenTruthAndRecoTrack, double &TempAngleBeetweenTruthAndRecoTrack)
int CountBadLeadingMuonTrackAndLeadingPlusSecondBadPurity
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
std::vector< TrackID > TrackIDs
const simb::MCParticle * TrackIdToParticle_P(int id) const
int GoodEvents4OrMoreMuonTrack
TH1D * h_Efficiency_ThetaXZ
TH2D * h_TrackTooShort_ThetaXZ_SinThetaYZ
TH2D * h_ThetaXZ_ThetaYZ_LeadingPlusSecondOk
int CountTrackLengthTooShort
TH2D * h_MuonTrackStitching_Distance_Angle
TH2D * h_MuonTrackStitching_FailedCriteria_CriteriaTwoTracks_Angle
Geometry information for a single TPC.
TH2D * h_NoMuonTrack_MaxTrackLength_PDGCode
TH2D * h_MuonTrackStitching_MatchedCriteria_TrackResLeading_TrackResSecond
TH2D * h_Efficiency_ThetaXZ_SinThetaYZ_LeadingPlusSecond
TH2D * h_Criteria_NRecoTrack_num
Vector_t VertexDirection() const
TH2D * h_MuonTrackStitching_FailedCriteria_Distance_Angle
bool get(SelectorBase const &, Handle< PROD > &result) const
TH2D * h_ThetaXZ_ThetaYZ_LeadingPlusSecondOk_num
std::set< const gar::rec::Track * > TrackList
TH2D * h_MuonTrackStitching_MatchedCriteria_TrackRes_Completeness
TH1D * h_Efficiency_ThetaYZ
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)
int CountBadLeadingMuonTrackAndOnlyOneMuonTrack
TH2D * h_ThetaXZ_SinThetaYZ_LeadingPlusSecondOk
TH2D * h_MuonTrackStitching_MatchedCriteria_CriteriaTwoTracks_Angle
TH2D * h_MuonTrackStitching_FailedCriteria_CompletenessSecondMuon_Angle
TH2D * h_MuonTrackStitching_FailedCriteria_CompletenessLeading_CompletenessSecond
bool insideFV(double vertex[4])
TH2D * h_FailedReconstruction_ThetaXZ_SinThetaYZ
geo::TPCID FindTPCAtPosition(double const worldLoc[3]) const
Returns the ID of the TPC at specified location.
TH2D * h_MuonTrackStitching_FailedCriteria_TrackRes_Completeness
std::string fMCTruthModuleLabel
bool getByLabel(std::string const &label, std::string const &instance, Handle< PROD > &result) const
TH2D * h_MuonTrackStitching_TrackResLeading_TrackResSecond
double Length(size_t p=0) const
Access to various track properties.
TH2D * h_MuonTrackStitching_FailedCriteriaAndLeadingPlusSecondBad_TrackResLeading_TrackResSecond
#define DEFINE_ART_MODULE(klass)
TH2D * h_Efficiency_ThetaXZ_ThetaYZ_DifferenceLeadingAndLeadingPlusSecond
T get(std::string const &key) const
double MCTruthMuonVertex[4]
int CountBadLeadingMuonTrackAndLeadingPlusSecondBadTrackTooShort
TH2D * h_Efficiency_ThetaXZ_ThetaYZ
TH2D * h_MuonTrackStitching_FailedCriteriaAndLeadingPlusSecondBad_TrackRes_Completeness
Point_t const & Vertex() const
int CountBadLeadingMuonTrack
void processEff(const art::Event &evt, bool &isFiducial)
int CountGoodLeadingMuonTrack
int CountBadLeadingMuonTrackAndOnlyOneMuonTrackTrackTooLong
TH2D * h_Criteria_NMuonTrack_den
void FuncDistanceAndAngleBetweenTracks(art::Ptr< recob::Track > Track1, art::Ptr< recob::Track > Track2, double &TempDistanceBetweenTracks, double &TempAngleBetweenTracks, double &TempCriteriaTwoTracks)
The data type to uniquely identify a TPC.
TH1D * h_Efficiency_SinThetaYZ
Definition of data types for geometry description.
TH2D * h_ThetaXZ_SinThetaYZ_den
Detector simulation of raw signals on wires.
TH2D * h_MuonTrackStitching_FailedCriteria_TrackResSecondMuon_Angle
const sim::ParticleList & ParticleList() const
double MCTruthMuonThetaYZ
double Vx(const int i=0) const
int CountBadLeadingMuonTrackAndLeadingPlusSecondBad
TH2D * h_Criteria_NMuonTrack
TH2D * h_Efficiency_ThetaXZ_ThetaYZ_LeadingPlusSecond
TH2D * h_MuonTrackStitching_FailedCriteriaButLeadingPlusSecondGood_TrackResLeading_TrackResSecond
TH2D * h_ThetaXZ_SinThetaYZ_LeadingPlusSecondOk_num
int CountBadLeadingMuonTrackAndLeadingPlusSecondBadTrackTooLong
Contains all timing reference information for the detector.
TH2D * h_Criteria_NRecoTrack_den
TH2D * h_Purity_ThetaXZ_SinThetaYZ
TH2D * h_ThetaXZ_ThetaYZ_num
TH2D * h_MuonTrackStitching_FailedCriteriaAndLeadingPlusSecondBad_Distance_Angle
Vector_t EndDirection() const
int CountTrackLengthTooLong
TH2D * h_ThetaXZ_ThetaYZ_den
const TLorentzVector & Momentum(const int i=0) const
TH2D * h_MuonTrackStitching_CompletenessSecondMuon_Angle
Provides recob::Track data product.
void analyze(const art::Event &evt) override
double Vz(const int i=0) const
TH2D * h_Efficiency_ThetaXZ_SinThetaYZ
Point_t const & End() const
TH2D * h_Criteria_NRecoTrack
TrackCollectionProxyElement< TrackCollProxy > Track
Proxy to an element of a proxy collection of recob::Track objects.
TH2D * h_MuonTrackStitching_FailedCriteriaButLeadingPlusSecondGood_TrackRes_Completeness
art::ServiceHandle< geo::Geometry const > geom
TH2D * h_FailedReconstruction_ThetaXZ_ThetaYZ
double truthLength(const simb::MCParticle *MCparticle)
TH2D * h_MuonTrackStitching_FailedCriteria_TrackResLeading_TrackResSecond
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
int CountBadLeadingMuonTrackAndOnlyOneMuonTrackTrackTooShort
TH2D * h_MuonTrackStitching_CriteriaTwoTracks_Angle
double MCTruthMuonThetaXZ
LArSoft geometry interface.
TH2D * h_Criteria_NMuonTrack_num
TH2D * h_MuonTrackStitching_TrackResSecondMuon_Angle
int CountBadLeadingMuonTrackAndOnlyOneMuonTrackCompleteness
int CountBadLeadingMuonTrackAndLeadingPlusSecondBadCompleteness
void LocalToWorld(const double *tpc, double *world) const
Transform point from local TPC frame to world frame.
double MCTruthMuonMomentum
TH2D * h_ThetaXZ_SinThetaYZ_num
double Vy(const int i=0) const
TH2D * h_MuonTrackStitching_MatchedCriteria_Distance_Angle
QTextStream & endl(QTextStream &s)
Event finding and building.
TH2D * h_MuonTrackStitching_FailedCriteriaButLeadingPlusSecondGood_Distance_Angle
std::string fTrackModuleLabel