1121 simb::MCParticle* MCTruthMuonParticle = NULL;
1123 art::ServiceHandle<cheat::ParticleInventoryService const> pi_serv;
1124 const sim::ParticleList& plist = pi_serv->ParticleList();
1125 simb::MCParticle* particle = 0;
1127 for (sim::ParticleList::const_iterator ipar = plist.begin(); ipar != plist.end(); ++ipar) {
1128 particle = ipar->second;
1130 if (particle->Mother() ==
1132 const TLorentzVector& positionStart = particle->Position(0);
1133 positionStart.GetXYZT(
1138 MCTruthMuonParticle = particle;
1141 sqrt(pow(particle->Momentum().Px(), 2) + pow(particle->Momentum().Py(), 2) +
1142 pow(particle->Momentum().Pz(), 2));
1144 if (particle->Momentum().Pz() >= 0 && particle->Momentum().Px() >= 0) {
1146 (180.0 / 3.14159) * atan(particle->Momentum().Px() / particle->Momentum().Pz());
1148 else if (particle->Momentum().Pz() < 0 && particle->Momentum().Px() >= 0) {
1150 180.0 + (180.0 / 3.14159) * atan(particle->Momentum().Px() / particle->Momentum().Pz());
1152 else if (particle->Momentum().Pz() < 0 && particle->Momentum().Px() < 0) {
1154 180.0 + (180.0 / 3.14159) * atan(particle->Momentum().Px() / particle->Momentum().Pz());
1156 else if (particle->Momentum().Pz() >= 0 && particle->Momentum().Px() < 0) {
1158 360.0 + (180.0 / 3.14159) * atan(particle->Momentum().Px() / particle->Momentum().Pz());
1165 double MCTruthLengthMuon =
truthLength(MCTruthMuonParticle);
1172 if (!isFiducial)
return;
1175 if (MCTruthMuonParticle) {
1189 int NMuonTracks = 0;
1191 art::Handle<std::vector<recob::Track>> TrackListHandle;
1193 std::vector<art::Ptr<recob::Track>> TrackList;
1194 art::fill_ptr_vector(TrackList, TrackListHandle);
1195 int NRecoTracks = TrackList.size();
1196 art::FindManyP<recob::Hit> track_hits(TrackListHandle, event,
fTrackModuleLabel);
1197 if (NRecoTracks == 0) {
1198 MF_LOG_DEBUG(
"MuonTrackingEff") <<
"There are no reco tracks... bye";
1199 std::cout <<
"There are no reco tracks! MCTruthMuonThetaXZ: " << std::endl;
1212 MF_LOG_DEBUG(
"MuonTrackingEff") <<
"Found this many reco tracks " << NRecoTracks;
1216 double PurityLeadingMuon = 0.;
1217 double CompletenessLeadingMuon = 0.;
1218 double RecoLengthLeadingMuon = 0.;
1219 art::Ptr<recob::Track> TrackLeadingMuon;
1221 double RecoLengthSecondMuon = 0.;
1222 double CompletenessSecondMuon = 0.;
1223 double PuritySecondMuon = 0.;
1224 art::Ptr<recob::Track> TrackSecondMuon;
1226 double TrackLengthMuonSum = 0.;
1227 double tmpTotalRecoEnergy = 0.;
1229 double MaxLengthNoRecoMuon = 0;
1230 int PDGCodeMaxLengthNoRecoMuon = 0;
1232 const simb::MCParticle* RecoMuonParticle = NULL;
1234 std::vector<art::Ptr<recob::Hit>> tmp_TrackHits = track_hits.at(0);
1235 std::vector<art::Ptr<recob::Hit>> AllHits;
1236 art::Handle<std::vector<recob::Hit>> HitHandle;
1237 if (event.get(tmp_TrackHits[0].id(), HitHandle)) art::fill_ptr_vector(AllHits, HitHandle);
1239 auto const clockData =
1240 art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(event);
1243 for (
int i = 0; i < NRecoTracks; i++) {
1244 art::Ptr<recob::Track>
track = TrackList[i];
1245 std::vector<art::Ptr<recob::Hit>> TrackHits = track_hits.at(i);
1246 double tmpPurity = 0.;
1247 double tmpCompleteness = 0.;
1248 const simb::MCParticle* particle;
1251 clockData, AllHits, TrackHits, particle, tmpPurity, tmpCompleteness, tmpTotalRecoEnergy);
1254 std::cout <<
"ERROR: Truth matcher didn't find a particle!" << std::endl;
1258 if (track->Length() > MaxLengthNoRecoMuon && particle->PdgCode() !=
fMuonPDGCode &&
1260 MaxLengthNoRecoMuon = track->Length();
1261 PDGCodeMaxLengthNoRecoMuon = particle->PdgCode();
1265 tmpCompleteness > 0 && tmpPurity > 0) {
1268 TrackLengthMuonSum += track->Length();
1270 if (NMuonTracks == 1) {
1271 CompletenessLeadingMuon = tmpCompleteness;
1272 PurityLeadingMuon = tmpPurity;
1273 RecoLengthLeadingMuon = track->Length();
1274 TrackLeadingMuon =
track;
1276 RecoMuonParticle = particle;
1279 if (NMuonTracks >= 2 && tmpCompleteness > CompletenessLeadingMuon) {
1281 CompletenessSecondMuon = CompletenessLeadingMuon;
1282 PuritySecondMuon = PurityLeadingMuon;
1283 RecoLengthSecondMuon = RecoLengthLeadingMuon;
1284 TrackSecondMuon = TrackLeadingMuon;
1286 CompletenessLeadingMuon = tmpCompleteness;
1287 PurityLeadingMuon = tmpPurity;
1288 RecoLengthLeadingMuon = track->Length();
1289 TrackLeadingMuon =
track;
1291 RecoMuonParticle = particle;
1294 else if (NMuonTracks >= 2 && tmpCompleteness < CompletenessLeadingMuon &&
1295 tmpCompleteness > CompletenessSecondMuon) {
1296 CompletenessSecondMuon = tmpCompleteness;
1297 PuritySecondMuon = tmpPurity;
1298 RecoLengthSecondMuon = track->Length();
1299 TrackSecondMuon =
track;
1309 if (RecoMuonParticle && MCTruthMuonParticle) {
1313 h_TrackRes->Fill(RecoLengthLeadingMuon / MCTruthLengthMuon);
1315 std::cout <<
"TrackLeadingMuon->Vertex().X(): " << TrackLeadingMuon->Vertex().X()
1317 std::cout <<
"MCTruthMuonParticle->Vz(): " << MCTruthMuonParticle->Vz() << std::endl;
1320 double DistanceBetweenTruthAndRecoTrack;
1321 double AngleBetweenTruthAndRecoTrack;
1324 DistanceBetweenTruthAndRecoTrack,
1325 AngleBetweenTruthAndRecoTrack);
1327 h_VertexRes->Fill(DistanceBetweenTruthAndRecoTrack);
1331 CompletenessLeadingMuon);
1332 if (NMuonTracks >= 2) {
1333 double DistanceBetweenTracks;
1334 double AngleBetweenTracks;
1335 double CriteriaTwoTracks;
1339 DistanceBetweenTracks,
1344 RecoLengthLeadingMuon / MCTruthLengthMuon, RecoLengthSecondMuon / MCTruthLengthMuon);
1347 RecoLengthSecondMuon / MCTruthLengthMuon, AngleBetweenTracks);
1349 AngleBetweenTracks);
1354 if (CompletenessLeadingMuon < 0.5) {
1363 if (PurityLeadingMuon < 0.5) {
1372 if (RecoLengthLeadingMuon / MCTruthLengthMuon < 0.75) {
1381 if (RecoLengthLeadingMuon / MCTruthLengthMuon > 1.25) {
1409 if (CompletenessLeadingMuon < 0.5 || PurityLeadingMuon < 0.5 ||
1410 RecoLengthLeadingMuon / MCTruthLengthMuon < 0.75 ||
1411 RecoLengthLeadingMuon / MCTruthLengthMuon > 1.25) {
1417 RecoLengthLeadingMuon / MCTruthLengthMuon, CompletenessLeadingMuon);
1424 if (NMuonTracks >= 2) {
1425 double AngleBetweenTracks;
1426 double DistanceBetweenTracks;
1427 double CriteriaTwoTracks;
1430 DistanceBetweenTracks,
1434 if (AngleBetweenTracks > 160.) {
1436 std::cout <<
"Angle b/w tracks: " << AngleBetweenTracks << std::endl;
1437 std::cout <<
"CriteriaTwoTracks: " << CriteriaTwoTracks << std::endl;
1438 std::cout <<
"CompletenessLeadingMuon: " << CompletenessLeadingMuon << std::endl;
1439 std::cout <<
"CompletenessSecondMuon: " << CompletenessSecondMuon << std::endl;
1440 std::cout <<
"PurityLeadingMuon: " << PurityLeadingMuon << std::endl;
1441 std::cout <<
"PuritySecondMuon: " << PuritySecondMuon << std::endl;
1442 std::cout <<
"TrackLeadingMuon: " << RecoLengthLeadingMuon / MCTruthLengthMuon
1444 std::cout <<
"TrackResSecondMuon: " << RecoLengthSecondMuon / MCTruthLengthMuon
1446 std::cout <<
"TrackID leading: " << TrackLeadingMuon->ID() << std::endl;
1447 std::cout <<
"TrackID second: " << TrackSecondMuon->ID() << std::endl;
1451 AngleBetweenTracks);
1453 RecoLengthLeadingMuon / MCTruthLengthMuon, RecoLengthSecondMuon / MCTruthLengthMuon);
1455 CompletenessLeadingMuon, CompletenessSecondMuon);
1457 RecoLengthSecondMuon / MCTruthLengthMuon, AngleBetweenTracks);
1459 CompletenessSecondMuon, AngleBetweenTracks);
1461 AngleBetweenTracks);
1462 if ((CompletenessLeadingMuon + CompletenessSecondMuon) >= 0.5 &&
1463 PurityLeadingMuon >= 0.5 && PuritySecondMuon >= 0.5 &&
1464 (RecoLengthLeadingMuon + RecoLengthSecondMuon) / MCTruthLengthMuon >= 0.75 &&
1465 (RecoLengthLeadingMuon + RecoLengthSecondMuon) / MCTruthLengthMuon <= 1.25) {
1472 RecoLengthLeadingMuon / MCTruthLengthMuon, CompletenessLeadingMuon);
1474 ->Fill(RecoLengthLeadingMuon / MCTruthLengthMuon,
1475 RecoLengthSecondMuon / MCTruthLengthMuon);
1477 DistanceBetweenTracks, AngleBetweenTracks);
1479 if ((CompletenessLeadingMuon + CompletenessSecondMuon) < 0.5 || PurityLeadingMuon < 0.5 ||
1480 PuritySecondMuon < 0.5 ||
1481 (RecoLengthLeadingMuon + RecoLengthSecondMuon) / MCTruthLengthMuon < 0.75 ||
1482 (RecoLengthLeadingMuon + RecoLengthSecondMuon) / MCTruthLengthMuon > 1.25) {
1486 RecoLengthLeadingMuon / MCTruthLengthMuon, CompletenessLeadingMuon);
1488 ->Fill(RecoLengthLeadingMuon / MCTruthLengthMuon,
1489 RecoLengthSecondMuon / MCTruthLengthMuon);
1491 DistanceBetweenTracks, AngleBetweenTracks);
1493 if ((CompletenessLeadingMuon + CompletenessSecondMuon) < 0.5) {
1496 if (PurityLeadingMuon < 0.5 || PuritySecondMuon < 0.5) {
1499 if ((RecoLengthLeadingMuon + RecoLengthSecondMuon) / MCTruthLengthMuon < 0.75) {
1502 if ((RecoLengthLeadingMuon + RecoLengthSecondMuon) / MCTruthLengthMuon > 1.25) {
1506 else if (NMuonTracks == 1) {
1508 if (CompletenessLeadingMuon < 0.5) {
1512 if (RecoLengthLeadingMuon / MCTruthLengthMuon < 0.75) {
1515 if (RecoLengthLeadingMuon / MCTruthLengthMuon > 1.25) {
1521 if (CompletenessLeadingMuon >= 0.5 && PurityLeadingMuon >= 0.5 &&
1522 RecoLengthLeadingMuon / MCTruthLengthMuon >= 0.75 &&
1523 RecoLengthLeadingMuon / MCTruthLengthMuon <= 1.25) {
1536 RecoLengthLeadingMuon / MCTruthLengthMuon, CompletenessLeadingMuon);
1543 if (NMuonTracks >= 2) {
1545 RecoLengthLeadingMuon / MCTruthLengthMuon, RecoLengthSecondMuon / MCTruthLengthMuon);
1547 double AngleBetweenTracks;
1548 double DistanceBetweenTracks;
1549 double CriteriaTwoTracks;
1552 DistanceBetweenTracks,
1556 AngleBetweenTracks);
1558 AngleBetweenTracks);
1563 if (!RecoMuonParticle && MCTruthMuonParticle) {
1575 static_cast<double>(PDGCodeMaxLengthNoRecoMuon));
TH2D * h_MuonTrackStitching_TrackRes_Completeness
TH2D * h_NoRecoTrackAtAll_ThetaXZ_SinThetaYZ
TH2D * h_NoMuonTrack_ThetaXZ_SinThetaYZ
TH2D * h_Completeness_ThetaXZ_SinThetaYZ
int CountBadLeadingMuonTrackAndOnlyOneMuonTrackPurity
int CountBadLeadingMuonTrackButLeadingPlusSecondGood
int BadEvents4OrMoreMuonTrack
TH2D * h_TrackTooLong_ThetaXZ_SinThetaYZ
void FuncDistanceAndAngleBetweenTruthAndRecoTrack(const simb::MCParticle *&MCparticle, art::Ptr< recob::Track > Track, double &TempDistanceBetweenTruthAndRecoTrack, double &TempAngleBeetweenTruthAndRecoTrack)
int CountBadLeadingMuonTrackAndLeadingPlusSecondBadPurity
int GoodEvents4OrMoreMuonTrack
TH2D * h_TrackTooShort_ThetaXZ_SinThetaYZ
TH2D * h_ThetaXZ_ThetaYZ_LeadingPlusSecondOk
int CountTrackLengthTooShort
TH2D * h_MuonTrackStitching_Distance_Angle
TH2D * h_MuonTrackStitching_FailedCriteria_CriteriaTwoTracks_Angle
TH2D * h_NoMuonTrack_MaxTrackLength_PDGCode
TH2D * h_MuonTrackStitching_MatchedCriteria_TrackResLeading_TrackResSecond
TH2D * h_Criteria_NRecoTrack_num
TH2D * h_MuonTrackStitching_FailedCriteria_Distance_Angle
TH2D * h_MuonTrackStitching_MatchedCriteria_TrackRes_Completeness
process_name use argoneut_mc_hitfinder track
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
TH2D * h_MuonTrackStitching_FailedCriteria_TrackRes_Completeness
TH2D * h_MuonTrackStitching_TrackResLeading_TrackResSecond
TH2D * h_MuonTrackStitching_FailedCriteriaAndLeadingPlusSecondBad_TrackResLeading_TrackResSecond
double MCTruthMuonVertex[4]
int CountBadLeadingMuonTrackAndLeadingPlusSecondBadTrackTooShort
TH2D * h_MuonTrackStitching_FailedCriteriaAndLeadingPlusSecondBad_TrackRes_Completeness
int CountBadLeadingMuonTrack
int CountGoodLeadingMuonTrack
int CountBadLeadingMuonTrackAndOnlyOneMuonTrackTrackTooLong
void FuncDistanceAndAngleBetweenTracks(art::Ptr< recob::Track > Track1, art::Ptr< recob::Track > Track2, double &TempDistanceBetweenTracks, double &TempAngleBetweenTracks, double &TempCriteriaTwoTracks)
TH2D * h_ThetaXZ_SinThetaYZ_den
TH2D * h_MuonTrackStitching_FailedCriteria_TrackResSecondMuon_Angle
double MCTruthMuonThetaYZ
int CountBadLeadingMuonTrackAndLeadingPlusSecondBad
TH2D * h_MuonTrackStitching_FailedCriteriaButLeadingPlusSecondGood_TrackResLeading_TrackResSecond
int CountBadLeadingMuonTrackAndLeadingPlusSecondBadTrackTooLong
TH2D * h_Purity_ThetaXZ_SinThetaYZ
TH2D * h_ThetaXZ_ThetaYZ_num
TH2D * h_MuonTrackStitching_FailedCriteriaAndLeadingPlusSecondBad_Distance_Angle
int CountTrackLengthTooLong
TH2D * h_ThetaXZ_ThetaYZ_den
TH2D * h_MuonTrackStitching_CompletenessSecondMuon_Angle
TH2D * h_MuonTrackStitching_FailedCriteriaButLeadingPlusSecondGood_TrackRes_Completeness
TH2D * h_FailedReconstruction_ThetaXZ_ThetaYZ
double truthLength(const simb::MCParticle *MCparticle)
TH2D * h_MuonTrackStitching_FailedCriteria_TrackResLeading_TrackResSecond
int CountBadLeadingMuonTrackAndOnlyOneMuonTrackTrackTooShort
TH2D * h_MuonTrackStitching_CriteriaTwoTracks_Angle
double MCTruthMuonThetaXZ
TH2D * h_Criteria_NMuonTrack_num
TH2D * h_MuonTrackStitching_TrackResSecondMuon_Angle
int CountBadLeadingMuonTrackAndOnlyOneMuonTrackCompleteness
int CountBadLeadingMuonTrackAndLeadingPlusSecondBadCompleteness
double MCTruthMuonMomentum
TH2D * h_ThetaXZ_SinThetaYZ_num
BEGIN_PROLOG could also be cout
TH2D * h_MuonTrackStitching_MatchedCriteria_Distance_Angle
TH2D * h_MuonTrackStitching_FailedCriteriaButLeadingPlusSecondGood_Distance_Angle
std::string fTrackModuleLabel