72 : fParams(params()), fHaveTruthInfo(false)
121 for(
size_t l = 0;
l < flslistHandle->size(); ++
l){
122 const std::vector<sim::FLSHit> &
hits = flslistHandle->at(
l).fHits;
123 for(
size_t f = 0;
f < hits.size(); ++
f){
146 for(
size_t i = 0;
i < phots->size(); ++
i){
158 mf::LogWarning(
"BackTracker") <<
"unable to find PhotonTransport " 160 <<
"any attempt to get the PhotonTransport objects from " 161 <<
"the backtracker will fail";
169 <<
"any attempt to get the Particles from " 170 <<
"the backtracker will fail";
174 std::set<art::Ptr<simb::MCTruth>> truths;
176 std::map<int, art::Ptr<simb::MCTruth>> trackIDToMCTruthPtr;
179 for(
size_t p = 0;
p < pHandle->size(); ++
p){
188 trackIDToMCTruthPtr.emplace(pHandle->at(
p).TrackId(), mct);
191 mf::LogWarning(
"BackTracker") <<
"unable to find MCTruth from Particle list " 193 <<
"any attempt to get the MCTruth objects from " 194 <<
"the backtracker will fail\n" 195 <<
"message from caught exception:\n" << ex;
203 std::map<art::Ptr<simb::MCTruth>,
int> truthIdxs;
209 for(
auto it: trackIDToMCTruthPtr){
228 for(
size_t p = 0;
p < pHandle->size(); ++
p){
229 int parTrId = pHandle->at(
p).TrackId();
232 bool match_truep =
false;
233 for(
size_t te = 0; te < teAssn.at(
p).size(); ++te){
234 if (parTrId == teAssn.at(
p)[te]->TrackId()) {
243 for(
size_t par2 = 0; par2 < pHandle->size(); ++par2){
244 for(
size_t te = 0; te < teAssn.at(par2).size(); ++te){
245 if (parTrId == teAssn.at(par2)[te]->TrackId()) {
250 if (match_truep)
break;
252 if (match_truep)
break;
268 mf::LogWarning(
"BackTracker") <<
"unable to find sim::TrueEnergy " 270 <<
"any attempt to get the sim::TrueEnergy's from " 271 <<
"the backtracker will fail";
278 unsigned int const&
cell)
const 283 LOG_DEBUG(
"BackTracker") <<
"no collection of sim::PhotonSignal for " 284 <<
" plane = " << plane <<
" cell = " << cell
285 <<
" probably this was a noise hit, returning an" 287 const std::vector<sim::PhotonSignal> empty;
300 LOG_DEBUG(
"BackTracker") <<
"no collection of sim::PhotonSignal for " 301 <<
" plane = " << hit->
Plane() <<
" cell = " << hit->
Cell()
302 <<
" probably this was a noise hit, returning an" 304 const std::vector<sim::PhotonSignal> empty;
309 const std::vector<sim::PhotonSignal>& phots =
it->second;
311 std::vector<sim::PhotonSignal> hitPhots;
315 for(
unsigned int i = 0;
i < phots.size(); ++
i){
318 hitPhots.push_back(phot);
326 unsigned int const&
cell)
const 331 LOG_DEBUG(
"BackTracker") <<
"no collection of sim::FLSHit for " 332 <<
" plane = " << plane <<
" cell = " << cell
333 <<
" probably this was a noise hit, returning an" 335 const std::vector<sim::FLSHit> empty;
338 return std::vector<sim::FLSHit>(
it->second.begin(),
it->second.end());
344 std::map<int, int> idToNPhot;
349 std::vector<std::pair<int, int> > nphotToId;
350 nphotToId.reserve(idToNPhot.size());
351 for(std::map<int, int>::iterator
it = idToNPhot.begin();
it != idToNPhot.end(); ++
it){
355 if(nphotToId.empty()){
356 LOG_DEBUG(
"BackTracker") <<
"No matching photons, this is probably a noise hit. Returning empty vector";
357 return std::vector<sim::FLSHit>();
362 std::sort(nphotToId.rbegin(), nphotToId.rend());
367 <<
"PhotonSignal with no matching FLSHit. That shouldn't be possible\n" 368 << __FILE__ <<
":" << __LINE__ <<
"\n";
371 std::vector<sim::FLSHit>
ret;
372 const std::list<sim::FLSHit>& flshits = it->second;
374 ret.reserve(nphotToId.size()*flshits.size());
376 for(
unsigned int idIdx = 0; idIdx < nphotToId.size(); ++idIdx){
377 const int trackid = nphotToId[idIdx].second;
379 if(fls.GetTrackID() == trackid){
397 const std::list<sim::FLSHit>& flshits =
it->second;
402 if(fls.GetTrackID() == trackId){
403 energy += ((useBirksE) ? fls.GetEdepBirks() : fls.GetEdep());
422 std::map<int, int> idToNPhot;
428 for(std::map<int, int>::iterator
it = idToNPhot.begin();
it != idToNPhot.end(); ++
it){
429 if(
it->second > mostPhot){
430 mostPhot =
it->second;
437 mf::LogWarning(
"BackTracker") <<
"No photons in this cell are in time with the hit. " 438 <<
"Returning null pointer";
464 for(
size_t h = 0;
h <
hits.size(); ++
h){
466 if (
IsNoise(hit)) totalNoise += 1;
469 int totalPhysics = total - totalNoise;
470 double physicsFrac = (double)totalPhysics / (
double)
total;
486 mf::LogWarning(
"BackTracker") <<
"can't find particle with track id " 487 <<
id <<
" in sim::ParticleNavigator" 488 <<
" returning null pointer";
492 return part_it->second;
512 <<
"attempting to find MCTruth index for " 513 <<
"out of range value: " << mct
515 << __FILE__ <<
":" << __LINE__ <<
"\n";
529 std::vector<const sim::Particle*>
ret;
531 for(std::set<int>::const_iterator itr =
fTrackIDs.begin(); itr !=
fTrackIDs.end(); itr++){
544 std::map<int, int> idToNPhot;
545 for(
size_t h = 0;
h < hits.size(); ++
h){
552 std::vector<std::pair<int, int> > nphotToId;
553 for(std::map<int, int>::iterator
it = idToNPhot.begin();
it != idToNPhot.end(); ++
it){
559 std::sort(nphotToId.rbegin(), nphotToId.rend());
561 std::vector<const sim::Particle*>
ret;
562 for(
unsigned int n = 0;
n < nphotToId.size(); ++
n){
564 if(part) ret.push_back(part);
577 std::map<int, double> sigtracks;
580 for(
size_t h = 0;
h < hits.size(); ++
h){
581 const std::vector<sim::FLSHit> fls =
HitToFLSHit(*hits[
h]);
582 for(
size_t f = 0;
f < fls.size(); ++
f){
584 sigtracks[fls[
f].GetTrackID()] += fls[
f].GetEdepBirks();
585 totalE += fls[
f].GetEdepBirks();
588 sigtracks[fls[
f].GetTrackID()] += fls[
f].GetEdep();
589 totalE += fls[
f].GetEdep();
595 std::vector<TrackIDE> trackIDE;
596 for(std::map<int, double>::iterator
it = sigtracks.begin();
it != sigtracks.end(); ++
it){
601 trackIDE.push_back(ide);
621 return HitsToTrackIDE(std::vector<const rb::CellHit*>(1, &hit), useBirksE);
632 const std::vector<sim::FLSHit> fhits = this->
HitToFLSHit(hit);
637 double weight = ((useBirksE) ? fhit.GetEdepBirks() : fhit.GetEdep());
640 x += weight * 0.5*(fhit.GetExitX() + fhit.GetEntryX());
641 y += weight * 0.5*(fhit.GetExitY() + fhit.GetEntryY());
642 z += weight * 0.5*(fhit.GetExitZ() + fhit.GetEntryZ());
648 if(w < 1.
e-5)
return TVector3(-999, -999, -999);
650 return TVector3(x/w, y/w, z/w);
656 const std::vector<rb::WeightedHit>& whits,
657 std::map<int, double>* purMap,
658 std::map<int, int>* parents,
659 bool energyPur)
const 661 const std::vector<RawWeightedHit> rwh(whits.begin(), whits.end());
668 const std::vector<RawWeightedHit>& whits,
669 std::map<int, double>* purMap,
670 std::map<int, int>* parents,
671 bool energyPur)
const 676 double total = whits.size();
682 if(purMap) purMap->clear();
684 for(
size_t h = 0;
h < whits.size(); ++
h){
686 double totalEInCell = 0;
687 double desiredEInCell = 0;
690 double weight = whits[
h].weight;
701 for(
size_t t = 0;
t < tracks.size(); ++
t){
702 const int trackID = tracks[
t].trackID;
703 if(trackIDs.find(trackID) != trackIDs.end()){
704 tally.insert(trackID);
707 std::map<int, int>::iterator parentsIt = parents->find(trackID);
708 if(parentsIt != parents->end()){
709 tally.insert(parentsIt->second);
716 for(
size_t t = 0;
t < tracks.size(); ++
t){
717 totalEInCell += tracks[
t].energy *
weight;
718 if(trackIDs.find(tracks[
t].trackID) != trackIDs.end()){
719 desiredE += tracks[
t].energy;
720 desiredEInCell += tracks[
t].energy;
728 if( desiredEInCell >= totalEInCell )
729 totalE += desiredEInCell;
731 totalE += totalEInCell;
734 for(std::set<int>::iterator
it = tally.begin();
it != tally.end(); ++
it){
739 if(!tally.empty()) ++desired;
746 if((total > 0) && !energyPur) purity = desired/total;
747 if((totalE > 0) && energyPur) purity = desiredE/totalE;
752 if(total > 0 && purMap){
753 for(std::map<int, double>::iterator
it = purMap->begin();
it != purMap->end(); ++
it){
771 const std::vector<const rb::CellHit*>&
hits,
772 std::map<int, double>* purMap,
773 std::map<int, int>* parents,
774 bool energyPur)
const 777 std::vector<RawWeightedHit> wHits;
787 const std::vector<const rb::CellHit*>&
hits,
788 const std::vector<const rb::CellHit*>& allhits,
790 std::map<int, double>* effMap,
792 double* desiredEnergy,
799 std::vector<RawWeightedHit> wHits;
803 effMap, energyEff, desiredEnergy,
804 totalEnergy, desiredHits, totHits);
809 const std::vector<rb::WeightedHit>& whits,
810 const std::vector<const rb::CellHit*>& allhits,
812 std::map<int, double>* effMap,
814 double* desiredEnergy,
819 const std::vector<RawWeightedHit> rwh(whits.begin(), whits.end());
821 effMap, energyEff, desiredEnergy, totalEnergy, desiredHits, totHits);
827 const std::vector<RawWeightedHit>& whits,
828 const std::vector<const rb::CellHit*>& allhits,
830 std::map<int, double>* effMap,
832 double* desiredEnergy,
840 assert(allhits.size() >= whits.size());
844 for(
unsigned int n = 0;
n < whits.size(); ++
n){
849 std::vector<const rb::CellHit*> hits_copy;
850 for(
int i = 0;
i < (
int)whits.size();
i++)
851 hits_copy.push_back(whits[
i].hit);
854 std::sort(hits_copy.begin(), hits_copy.end());
855 for(
unsigned int n = 0;
n+1 < hits_copy.size(); ++
n)
856 assert(hits_copy[
n] < hits_copy[
n+1]);
864 double desiredE = 0.;
867 if(effMap) effMap->clear();
868 std::map<int, double> totMap;
870 for(
size_t h = 0;
h < whits.size(); ++
h){
872 double totalEInCell = 0.;
873 double desiredEInCell = 0.;
876 double weight = whits[
h].weight;
881 bool hitCounted =
false;
889 for(
size_t t = 0;
t < tracks.size(); ++
t){
890 totalEInCell += tracks[
t].energy;
891 if(trackIDs.find(tracks[
t].trackID) != trackIDs.end() ){
892 desiredEInCell += tracks[
t].energy;
898 if(effMap) ++(*effMap)[tracks[
t].trackID];
903 if( (weight*totalEInCell) >= desiredEInCell)
904 desiredE += desiredEInCell;
906 desiredE += (weight*totalEInCell);
911 for(
size_t h = 0;
h < allhits.size(); ++
h){
917 bool hitCounted =
false;
923 for(
size_t t = 0;
t < tracks.size(); ++
t){
928 if(trackIDs.find(tracks[
t].trackID) != trackIDs.end()){
929 totalE += tracks[
t].energy;
935 if(effMap) ++totMap[tracks[
t].trackID];
942 if((total > 0.) && !energyEff) efficiency = desired/total;
943 if((totalE > 0.) && energyEff) efficiency = desiredE/totalE;
945 if((total > 0.) && desiredHits ) *desiredHits = desired;
946 if((total > 0.) && totHits ) *totHits = total;
947 if((totalE > 0.) && desiredEnergy) *desiredEnergy = desiredE;
948 if((totalE > 0.) &&
totalEnergy ) *totalEnergy = totalE;
951 for(std::map<int, double>::iterator
it = effMap->begin();
it != effMap->end(); ++
it){
952 if(totMap[
it->first] > 0)
953 it->second /= totMap[
it->first];
971 for(
size_t h = 0;
h < hits.
size(); ++
h){
975 for(
size_t t = 0;
t < tracks.size(); ++
t){
976 if(tracks[
t].trackID == trackID){
991 if ((nxcells >= 3) && (nycells >= 3))
return 1;
1007 return a.second.efficiency>b.second.efficiency;
1011 return a.second.purity>b.second.purity;
1017 const std::vector<const rb::CellHit*>& allHits,
1020 std::vector < std::pair<unsigned int,NeutrinoEffPur> > sliceNeutrinoInteractions;
1028 std::set<int> neutrinoTrackIDs;
1029 for (
unsigned int i = 0;
i < particleList.size(); ++
i){
1030 neutrinoTrackIDs.insert(particleList[
i]->TrackId());
1033 double energySlice = -1.0;
1034 double energyTotal = -1.0;
1042 NeutrinoEffPur sliceNeutrinoInteraction = {neutrinoInteraction, sliceEff, slicePur,
s, energySlice, energyTotal, nSliceHits, nTotalHits};
1043 sliceNeutrinoInteractions.push_back(
std::make_pair(
s,sliceNeutrinoInteraction) );
1047 if (!sortPur) std::sort(sliceNeutrinoInteractions.begin(),sliceNeutrinoInteractions.end(),
sort_idxeff);
1048 else std::sort(sliceNeutrinoInteractions.begin(),sliceNeutrinoInteractions.end(),
sort_idxpur);
1050 std::vector<int> nuIndex;
1052 for (
unsigned int i = 0;
i<sliceNeutrinoInteractions.size(); ++
i ){
1053 nuIndex.push_back(
int(sliceNeutrinoInteractions[
i].first) );
1063 const std::vector<const rb::CellHit*>& allHits,
1066 std::vector <NeutrinoEffPur> sliceNeutrinoInteractions;
1074 std::set<int> neutrinoTrackIDs;
1075 for (
unsigned int i = 0;
i < particleList.size(); ++
i){
1076 neutrinoTrackIDs.insert(particleList[
i]->TrackId());
1079 double energySlice = -1.0;
1080 double energyTotal = -1.0;
1088 NeutrinoEffPur sliceNeutrinoInteraction = {neutrinoInteraction, sliceEff, slicePur,
s, energySlice, energyTotal, nSliceHits, nTotalHits};
1089 sliceNeutrinoInteractions.push_back(sliceNeutrinoInteraction);
1093 if (!sortPur) std::sort(sliceNeutrinoInteractions.begin(),sliceNeutrinoInteractions.end(),
sort_eff);
1094 else std::sort(sliceNeutrinoInteractions.begin(),sliceNeutrinoInteractions.end(),
sort_pur);
1096 return sliceNeutrinoInteractions;
1102 const std::vector<const rb::CellHit*>& allHits,
1112 const std::vector<std::vector<cheat::NeutrinoEffPur>>& slTruthTable,
1117 for(
auto slc: slTruthTable )
1118 assert(
slc.size() == nusWithIdx.size() &&
"Badly shaped Slice/truth table!" );
1121 std::vector<int> nuIds( slTruthTable.size(), -1 );
1122 if( nusWithIdx.empty() )
return nuIds;
1125 std::vector<int> bestSlices( nusWithIdx.size() );
1126 for(
size_t nuIdx = 0; nuIdx < nusWithIdx.size(); ++nuIdx ){
1130 for(
size_t sliceIdx = 0; sliceIdx < slTruthTable.size(); ++sliceIdx ){
1134 if( slMetric(effpur) > maxValue ){
1135 maxValue = slMetric(effpur);
1136 bestSlice = sliceIdx;
1142 bestSlices[nuIdx] = bestSlice;
1147 for(
size_t sliceIdx = 0; sliceIdx < slTruthTable.size(); ++sliceIdx ){
1151 for(
size_t nuIdx = 0; nuIdx < slTruthTable[sliceIdx].size(); ++nuIdx ){
1155 if( nuMetric(effpur) > maxValue ){
1156 maxValue = nuMetric(effpur);
1162 if( bestNu >= 0 && bestSlices[bestNu] == (
int)sliceIdx )
1163 nuIds[sliceIdx] = bestNu;
1173 const std::vector<std::vector<cheat::NeutrinoEffPur>>& slTruthTable)
const 1183 const std::vector<std::vector<cheat::NeutrinoEffPur>>& slTruthTable)
const 1193 const std::vector<std::vector<cheat::NeutrinoEffPur>>& slTruthTable)
const 1203 const std::vector<std::vector<cheat::NeutrinoEffPur>>& slTruthTable)
const 1213 const std::vector<const rb::CellHit*>& allHits,
1216 std::vector <NeutrinoEffPur> sliceNeutrinoInteractions;
1222 std::set<int> neutrinoTrackIDs;
1223 for (
unsigned int i = 0;
i < particleList.size(); ++
i){
1224 neutrinoTrackIDs.insert(particleList[
i]->TrackId());
1227 double energySlice = -1.0;
1228 double energyTotal = -1.0;
1236 NeutrinoEffPur sliceNeutrinoInteraction = {neutrinoInteraction, sliceEff, slicePur,
s, energySlice, energyTotal, nSliceHits, nTotalHits};
1237 sliceNeutrinoInteractions.push_back(sliceNeutrinoInteraction);
1241 if (!sortPur) std::sort(sliceNeutrinoInteractions.begin(),sliceNeutrinoInteractions.end(),
sort_eff);
1242 else std::sort(sliceNeutrinoInteractions.begin(),sliceNeutrinoInteractions.end(),
sort_pur);
1244 return sliceNeutrinoInteractions;
1250 const std::vector<const rb::CellHit*>& allHits,
1258 std::vector< std::vector<cheat::NeutrinoEffPur> >
BackTracker:: 1262 unsigned int numSlices = sliceList.size();
1263 std::vector<NeutrinoWithIndex> truthList =
allMCTruth();
1264 unsigned int numTruths = truthList.size();
1267 std::vector< std::vector<cheat::NeutrinoEffPur> >
table(numSlices,std::vector<cheat::NeutrinoEffPur>(numTruths));
1270 for (
unsigned int s = 0;
s < numSlices;
s++){
1271 for (
unsigned int t = 0;
t < numTruths;
t++){
1273 table[
s][
t] = truthFill;
1277 for (
unsigned int s = 0;
s < numSlices;
s++){
1279 std::vector<double> sliceEnergies(numTruths,0.0);
1281 std::vector<int> sliceHits(numTruths,0);
1283 for (
unsigned int h = 0;
h < sliceList[
s]->NCell();
h++){
1285 if (
IsNoise(sliceList[
s]->Cell(
h)))
continue;
1287 std::vector<double> hitEnergies(numTruths,0.0);
1289 const std::vector<TrackIDE> trackIDs =
HitToTrackIDE(sliceList[
s]->Cell(
h));
1290 for (
unsigned int id = 0;
id < trackIDs.size();
id++){
1294 hitEnergies[truthID] += trackIDs[
id].energy;
1298 for (
unsigned int q = 0; q < numTruths; q++){
1299 if (hitEnergies[q] > 0.0){
1300 sliceEnergies[q] += hitEnergies[q];
1307 double totalSliceE = std::accumulate(sliceEnergies.begin(),sliceEnergies.end(),0.0);
1308 for (
unsigned int t = 0;
t < numTruths;
t++){
1311 if (totalSliceE > 0){
1312 purity = (sliceEnergies[
t])/totalSliceE;
1315 if ((purity < 0) || (purity > 1)) {
1317 <<
"The purity value "<<purity<<
" is not allowed. \n";
1321 table[
s][
t].energySlice = sliceEnergies[
t];
1322 table[
s][
t].nSliceHits = sliceHits[
t];
1327 for (
unsigned int t = 0;
t < numTruths;
t++){
1329 double totalMCTruthE = 0.0;
1330 int totalMCTruthHits = 0;
1332 for (
unsigned int s = 0;
s < numSlices;
s++){
1333 totalMCTruthE += table[
s][
t].energySlice;
1334 totalMCTruthHits += table[
s][
t].nSliceHits;
1337 for (
unsigned int s = 0;
s < numSlices;
s++){
1340 if (totalMCTruthE > 0){
1341 efficiency = (table[
s][
t].energySlice)/totalMCTruthE;
1344 if ((efficiency < 0) || (efficiency > 1)) {
1346 <<
"The efficiency value "<<efficiency<<
" is not allowed. \n";
1350 table[
s][
t].energyTotal = totalMCTruthE;
1351 table[
s][
t].nTotalHits = totalMCTruthHits;
1362 const std::set<int>& trackIDs,
1363 const std::set<int>& allowedDaughterIDs,
1364 const std::vector<const rb::CellHit*>& allHits)
const 1366 std::vector<ParticleEffPur>
matches;
1368 std::set<int> trkIDs = trackIDs;
1369 std::set<int> allowedDaughters = allowedDaughterIDs;
1371 std::map<int, int> parents;
1372 for(std::set<int>::iterator
it = trkIDs.begin();
it != trkIDs.end(); ++
it){
1374 if(!mother)
continue;
1385 for(
unsigned int itrack = 0; itrack < tracks.size(); ++itrack){
1387 std::map<int, double> purs;
1391 double bestpur = 0.0;
1392 for(std::set<int>::iterator
it = trkIDs.begin();
it != trkIDs.end(); ++
it){
1393 const double p = purs[*
it];
1401 double recoPur = bestpur;
1402 int recotrkid = bestid;
1404 double recoEff = 0.0;
1408 std::set<int>
ideff;
1409 ideff.insert(recotrkid);
1414 trkIDs.erase(trkIDs.find(recotrkid));
1416 std::map<int, int>::iterator parentsIt = parents.begin();
1417 while(parentsIt != parents.end()){
1418 if(parentsIt->second == recotrkid){ parents.erase(parentsIt++); }
1419 else{ ++parentsIt; }
1424 matches.push_back(match);
1442 allNeutrinoInteractions.push_back(neutrino);
1451 std::vector <NeutrinoWithIndex> allTruth;
1456 allTruth.push_back(neutrino);
1472 double kWindowNs = 1000;
1475 double kOffsetNs = 350;
1489 return fabs(dt+kOffsetNs) <= kWindowNs;
1494 std::map<int, int> & idToNPhot)
const 1501 const std::vector<sim::PhotonSignal>& phots =
it->second;
1505 for(
unsigned int i = 0;
i < phots.size(); ++
i){
1515 std::unordered_map< int, std::list<sim::FLSHit> >::const_iterator
it =
fTrackIDToFLSHit.find(trackID);
1517 std::vector<sim::FLSHit> empty;
1518 LOG_DEBUG(
"BackTracker") <<
"no collection of sim::FLSHit for " 1519 <<
"Track ID "<<trackID
1520 <<
". Returning empty vector";
1525 return std::vector<sim::FLSHit>(it->second.begin(), it->second.end());
1532 const std::vector<sim::FLSHit>& flshits_with_same_track_id =
ParticleToFLSHit(trackID);
1534 std::vector<sim::FLSHit> out_flshits;
1536 const int nflshits = flshits_with_same_track_id.size();
1539 for(
int i=0;
i<nflshits; ++
i){
1540 const sim::FLSHit& flshit= flshits_with_same_track_id[
i];
1542 out_flshits.push_back(flshit);
1552 std::vector<rb::Cluster> clusters(
fMCTruthList.size()+1);
1554 for(
unsigned int hiti = 0; hiti <
hits.size(); ++hiti) {
1565 clusters[
idx].Add(hit);
1582 unsigned int MCTidx;
1583 for(MCTidx = 0; MCTidx <
fMCTruthList.size(); ++MCTidx) {
1588 for(
unsigned int hiti = 0; hiti <
hits.size(); ++hiti) {
1592 if(this->
IsNoise(hit))
continue;
1596 if(idx == MCTidx) cluster.
Add(hit);
1606 const int& sliceIdx)
1611 unsigned int nsli = sliceList.size();
1613 for (
unsigned int isli = 0; isli < nsli; isli++){
1614 unsigned int nh = sliceList[isli]->NCell();
1615 for (
unsigned int ih = 0; ih < nh; ih++){
1625 for(
auto ide : ides)
1636 const std::vector<const rb::CellHit*>&
hits,
1637 const int& sliceIdx)
1644 double totalE = 0.0;
1645 unsigned int nh = hits.size();
1651 std::map< int, double> idToEnergyInProng;
1652 for(
unsigned int ih = 0; ih < nh; ih++){
1656 for(
auto ide: ides){
1657 idToEnergyInProng[ide.trackID] += ide.energy;
1658 totalE += ide.energy;
1668 if(idToEnergyInProng.size() == 0)
1671 for(
auto idtoe: idToEnergyInProng){
1672 if(idtoe.second > maxE ){
1673 maxE = idtoe.second;
1674 maxID = idtoe.first;
1679 puff.
purity = maxE/totalE;
1708 EnEsc = EnEnt = EnDep = -5;
1713 int TrajEnter, TrajExit;
1714 TrajEnter = TrajExit = -1;
1724 EnEnt = (
float)Par.
E( TrajEnter ) - Par.
Mass();
1728 EnEsc = (
float)Par.
E( TrajExit ) - Par.
Mass();
1731 EnDep = EnEnt - EnEsc;
1747 mf::LogWarning(
"BackTracker") <<
"Didn't find TrackId " << Par.
TrackId() <<
" in my map of sim::TrueEnergy's, returning -5.";
1758 for (
int dau=0; dau<NumDaught; ++dau) {
1761 DauEsc +=
std::max((
float)0., ThisEsc);
1771 DauEsc += DaughtEsc;
1794 float TotEnEsc =
std::max((
float)0., ParEEsc) +
std::max((
float)0., DauEEsc);
1795 TotEnDep = ParEDep - DauEEsc;
1797 if ( TotEnDep < 0 ) {
1804 if ( ParEDep != 0 &&
1816 TotEnDep =
std::max( 0., ParEDep - DauEEsc + Par.
Mass() );
1838 mf::LogWarning(
"BackTracker") <<
"Didn't find TrackId " << Par.
TrackId() <<
" in my map of sim::TrueEnergy's, returning -5.";
1855 TrajEnt = TrajEx = -5;
1860 bool EnteredVol =
false;
1861 bool PrevHit =
false;
1862 for (
unsigned int trj=0; trj<NumTraj; ++trj) {
1867 if ( !EnteredVol ) {
1872 if ( EnteredVol && InDet && (trj==NumTraj-1) ) {
1877 if ( PrevHit && !InDet) {
double E(const int i=0) const
::xsd::cxx::tree::id< char, ncname > id
T max(const caf::Proxy< T > &a, T b)
#define LOG_DEBUG(stream)
unsigned int NumberTrajectoryPoints() const
back track the reconstruction to the simulation
std::vector< const T * > PtrVecToVecRawPtr(const art::PtrVector< T > &xs) const
Helper function for implementing overloads.
std::unordered_map< geo::OfflineChan, std::vector< sim::PhotonSignal > > fCellToPhotonSignal
map of PhotonSignals in each cell in event
int GetPlaneID() const
Plane ID.
int32_t TDC() const
The time of the last baseline sample.
BackTracker(const Parameters ¶ms, art::ActivityRegistry ®)
fvar< T > fabs(const fvar< T > &x)
std::set< int > fTrackIDs
G4 track ids for all particles in the events.
const sim::Particle * TrackIDToMotherParticle(int const &id) const
std::vector< TrackIDE > HitToTrackIDE(const rb::CellHit &hit, bool useBirksE=false) const
Convenience function. HitsToTrackIDE but for a single hit.
static bool sort_idxeff(const std::pair< int, NeutrinoEffPur > &a, const std::pair< int, NeutrinoEffPur > &b)
Tool for function SliceToNeutrinoIndex, tells how to sort by putting structures with highest efficien...
Atom< int > MinPhysicsHits
std::vector< NeutrinoEffPur > SliceToNeutrinoInteractions(const std::vector< const rb::CellHit * > &sliceHits, const std::vector< const rb::CellHit * > &allHits, bool sortPur=false) const
Given a collection of hits (often a slice), returns vector of structures of neutrino interactions...
void HitToParticlesMap(const std::vector< const rb::Cluster * > &sliceList, const int &sliceIdx)
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
Atom< std::string > GeantModuleLabel
unsigned short Plane() const
#define DEFINE_ART_SERVICE(svc)
float energyFrac
fraction of hit energy from the particle with this trackID
bool InterceptsDetector(const sim::Particle &Par, int &TrajEnt, int &TrajEx)
A function which returns the trajectory point at which the particle first enters and leave the detect...
int GetCellID() const
Cell ID.
std::vector< int > SliceToOrderedNuIdsByEff(const std::vector< cheat::NeutrinoWithIndex > &nusWithIdx, const std::vector< std::vector< cheat::NeutrinoEffPur >> &slTruthTable) const
std::vector< int > SliceToOrderedNuIdsByEffPur(const std::vector< cheat::NeutrinoWithIndex > &nusWithIdx, const std::vector< std::vector< cheat::NeutrinoEffPur >> &slTruthTable) const
std::vector< NeutrinoWithIndex > allMCTruth() const
const double EnergyFromTrackId(const art::Ptr< rb::CellHit > &hit, int trackId, bool useBirksE=false) const
The total energy deposited in this cell by all FLSHits with this trackId.
std::vector< sim::FLSHit > ParticleToFLSHit(const int &trackID) const
All the FLSHits that were created by the track id trackID, sorted from most to least light...
std::map< double, std::vector< cheat::TrackIDE > > fCellToIDs
cheat::ParticleEffPur ClusterToParticle(const std::vector< const rb::Cluster * > &sliceList, const std::vector< const rb::CellHit * > &hits, const int &sliceIdx)
list_type::const_iterator const_iterator
Vertical planes which measure X.
int trackID
Truth particle track ID.
float energy
energy from the particle with this trackID
Atom< std::string > PhotonModuleLabel
static bool sort_idxpur(const std::pair< int, NeutrinoEffPur > &a, const std::pair< int, NeutrinoEffPur > &b)
Tool for function SliceToNeutrinoIndex, tells how to sort by putting structures with highest purity f...
A collection of associated CellHits.
std::vector< NeutrinoWithIndex > allNeutrinoInteractions() const
Function primarily for CAFMaker - returns a vector of all MCTruth interactions along with their truth...
std::vector< rb::Cluster > ClusterByTruth(std::vector< art::Ptr< rb::CellHit > > const &hits)
:
float CalcEscapingEnergy(const sim::Particle &Par, bool UseMCPart=true, std::string G4Label="trueens")
A function to calculate the amount of energy from a given MCParticle which is not contained in the de...
::xsd::cxx::tree::exception< char > exception
int fNumTrueEnergyWarnings
float CalcDaughterEscapingEnergy(const sim::Particle &Par)
A function to be used in conjunction with CalcTotalEscapingEnergy to calculate the energy which the d...
A single unit of energy deposition in the liquid scintillator.
std::pair< Spectrum *, CheatDecomp * > make_pair(SpectrumLoaderBase &loader_data, SpectrumLoaderBase &loader_mc, HistAxis *axis, Cut *cut, const SystShifts &shift, const Var &wei)
double purity
Putity of particle relative to track.
const sim::Particle * HitToParticle(art::Ptr< rb::CellHit > const &hit, bool quiet=false) const
Returns the sim::Particle object contributing the most light to the specified rb::CellHit.
std::vector< int > SliceToOrderedNuIdsByPur(const std::vector< cheat::NeutrinoWithIndex > &nusWithIdx, const std::vector< std::vector< cheat::NeutrinoEffPur >> &slTruthTable) const
std::vector< std::vector< cheat::NeutrinoEffPur > > SlicesToMCTruthsTable(const std::vector< const rb::Cluster * > &sliceList) const
Given ALL the slices in an event, including the noise slice, returns a vector of vector of structures...
void AccumulateHitContributions(rb::CellHit const &hit, std::map< int, int > &idToNPhot) const
Internal helper, adds the contributions from each particle in this cell to the map.
double efficiency
Efficiency (based on FLS energy) of neutrino interaction relative to slice.
int NumberDaughters() const
Horizontal planes which measure Y.
int maxID(double arrayInput[])
float TotalEscEnergy() const
art::PtrVector< rb::CellHit > AllCells() const
Get all cells from both views.
bool IsNoise(const art::Ptr< rb::CellHit > &hit) const
Is this hit not associated with any particles?
double HitCollectionEfficiency(const std::set< int > &trackIDs, const std::vector< rb::WeightedHit > &whits, const std::vector< const rb::CellHit * > &allhits, const geo::View_t &view, std::map< int, double > *effMap=0, bool energyEff=false, double *desiredEnergy=0, double *totalEnergy=0, int *desiredHits=0, int *totalHits=0) const
Returns the fraction of all energy in an event from a specific set of Geant4 track IDs that are repre...
std::vector< int > SliceToOrderedNuIds(const std::vector< cheat::NeutrinoWithIndex > &nusWithIdx, const std::vector< std::vector< cheat::NeutrinoEffPur >> &slTruthTable, std::function< double(const cheat::NeutrinoEffPur &)> slMetric, std::function< double(const cheat::NeutrinoEffPur &)> nuMetric) const
Given a vector of indexed neutrino interaction and a vector of slice truth tables, returns a vector of neutrino interaction indices ordered by best match to the corresponding slice. Here, best match is determined according to the given cheat::NeutrinoEffPur functions for the slice and the nu.
int Daughter(const int i) const
std::unordered_map< geo::OfflineChan, std::list< sim::FLSHit > > fCellToFLSHit
map of FLSHits in each cell in event
Atom< std::string > TrueEnergyModuleLabel
virtual void Add(const art::Ptr< rb::CellHit > &cell, double weight=1)
std::vector< art::Ptr< simb::MCTruth > > fMCTruthList
all the MCTruths for this event
std::vector< int > SliceToOrderedNuIdsByEnergy(const std::vector< cheat::NeutrinoWithIndex > &nusWithIdx, const std::vector< std::vector< cheat::NeutrinoEffPur >> &slTruthTable) const
unsigned short Cell() const
std::unordered_map< int, std::list< sim::FLSHit > > fTrackIDToFLSHit
map of G4 track ids to FLSHits
std::map< int, float > fIdToEnergy
Far Detector at Ash River, MN.
int trackID
Geant4 supplied trackID.
bool IsHitsNoise(std::vector< art::Ptr< rb::CellHit > > const &hits) const
Of this collection of hits, is it mostly noise or physics hits? True indicates noise.
int PdgCode(Resonance_t res, int Q)
(resonance id, charge) -> PDG code
std::vector< art::Ptr< sim::TrueEnergy > > fTrueEnergyList
all the TrueEnergy's for this event
sim::ParticleNavigator fParticleNav
Particle navigator to map track ID to Particle.
std::vector< sim::FLSHit > HitToFLSHit(const rb::CellHit &hit) const
All the FLSHits that contributed to this hit, sorted from most to least light.
novadaq::cnv::DetId DetId() const
Prefer ds::DetectorService::DetId() instead.
double efficiency
Efficiency of particle relative to track.
Atom< bool > MCOverlayMode
bool CompareByEnergy(const TrackIDE &a, const TrackIDE &b)
Does a have less energy than b?
Near Detector in the NuMI cavern.
const art::Ptr< simb::MCTruth > & ParticleToMCTruth(const sim::Particle *p) const
double EnergyMetric(const cheat::NeutrinoEffPur &ep)
Function for NeutrinoEffPur's nu interaction to slice energy.
bool isInsideDetectorBigBox(const double x_cm, const double y_cm, const double z_cm) const
Is the particle inside the detector Big Box?
int fSliceIDForHitParticleMap
std::map< int, int > fTrackIDToTrueEIndex
map of track id values to TrueEnergy entry
static bool sort_pur(const NeutrinoEffPur &a, const NeutrinoEffPur &b)
Tool for function SliceToNeutrinoInteractions, tells how to sort by putting structures with highest p...
double HitCollectionPurity(const std::set< int > &trackIDs, const std::vector< rb::WeightedHit > &whits, std::map< int, double > *purMap=0, std::map< int, int > *parents=0, bool energyPur=false) const
Returns the fraction of hits in a collection that come from the specified Geant4 track ids...
const sim::Particle * TrackIDToParticle(int const &id) const
Returns a pointer to the sim::Particle object corresponding to the given TrackID. ...
bool CompareByNPhoton(const sim::PhotonSignal &a, const sim::PhotonSignal &b)
Helper for SortByNPhoton.
void Rebuild(const art::Event &evt)
std::vector< TrackIDE > HitsToTrackIDE(const std::vector< const rb::CellHit * > &hits, bool useBirksE=false) const
Returns vector of TrackIDE structs contributing to the given collection of hits.
double Vx(const int i=0) const
double EffPurMetric(const cheat::NeutrinoEffPur &ep)
Function for NeutrinoEffPur's nu interaction to slice efficiency * purity.
const std::vector< sim::PhotonSignal > HitToPhotonSignal(const art::Ptr< rb::CellHit > &hit) const
Returns the PhotonSignals contributing the signal in the specified hit.
const std::vector< sim::PhotonSignal > CellToPhotonSignal(unsigned int const &plane, unsigned int const &cell) const
Returns the PhotonSignals contributing the signal in the specified cell. WARNING: Use with extreme ca...
A rawdata::RawDigit with channel information decoded.
code to link reconstructed objects back to the MC truth information
std::vector< int > SliceToNeutrinoIndex(const std::vector< const rb::CellHit * > &sliceHits, const std::vector< const rb::CellHit * > &allHits, bool sortPur=false) const
Given a collection of hits (often a slice), returns vector of neutrino indices corresponding to the v...
int GetTrackID() const
Track ID.
const art::Ptr< simb::MCTruth > & TrackIDToMCTruth(int const &id) const
bool fHaveTruthInfo
Set by Rebuild.
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
GlobalSignal< detail::SignalResponseType::FIFO, void(Event const &)> sPreProcessEvent
std::map< int, int > fTrackIDToMCTruthIndex
map of track id values to MCTruthList entry
static bool sortFLSHit(sim::FLSHit const &a, sim::FLSHit const &b)
double energySlice
Sum of FLS hits from the neutrino contributing to hits included in the slice.
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
double Vz(const int i=0) const
assert(nhit_max >=nhit_nbins)
double PurMetric(const cheat::NeutrinoEffPur &ep)
Function for NeutrinoEffPur's nu interaction to slice purity.
static bool sort_eff(const NeutrinoEffPur &a, const NeutrinoEffPur &b)
Tool for function SliceToNeutrinoInteractions, tells how to sort by putting structures with highest e...
double EffMetric(const cheat::NeutrinoEffPur &ep)
Function for NeutrinoEffPur's nu interaction to slice efficiency.
const std::vector< sim::FLSHit > CellToFLSHit(unsigned int const &plane, unsigned int const &cell) const
Returns the FLSHits contributing the signal in the specified cell. WARNING: Use with extreme caution!...
bool PassMuonCuts(int trackID, art::PtrVector< rb::CellHit > const &hits) const
Tool for NumuEAna that allows one to see if primary muon (or any track ID you feed the function) cont...
static Var totalEnergy(const std::shared_ptr< CAFAnaModel > &model)
float GetEdep() const
Get total Energy deposited into the cell for the whole FLSHit.
ParticleNavigator Add(const int &offset) const
TVector3 HitToXYZ(art::Ptr< rb::CellHit > const &hit, bool useBirksE=false) const
Returns the XYZ position of the energy deposition for a given hit.
float TotalDepEnergy() const
std::vector< ParticleEffPur > TracksToParticles(const std::vector< const rb::Track * > &tracks, const std::set< int > &trkIDs, const std::set< int > &allowedDaughters, const std::vector< const rb::CellHit * > &allHits) const
Given a collection of reconstructed tracks, this method finds the best match for each reco track to t...
iterator find(const key_type &key)
novadaq::cnv::DetId fDetID
detector id
std::vector< NeutrinoEffPur > SliceToMCTruth(const std::vector< const rb::CellHit * > &sliceHits, const std::vector< const rb::CellHit * > &allHits, bool sortPur=false) const
Given a collection of hits (often a slice), returns vector of structures of MCTruth, efficiency, and purity of that neutrino interaction relative to the slice. Efficiency is defined as FLS energy from neutrino interaction in slice / total FLS energy from neutrino interaction in event. This vector is sorted from the highest efficiency interaction to lowest. This function returns all MCTruth, including those without neutrino interactions, i.e. cosmics.
Atom< double > MinContribFrac
float CalcTotalEscapingEnergy(const sim::Particle &Par, bool UseMCPart=true, std::string G4Label="trueens")
An extension of CalcEscapingEnergy which also takes into account the energy which the daughters of gi...
rb::Cluster MCTruthToCluster(art::Ptr< simb::MCTruth > const &truth, std::vector< art::Ptr< rb::CellHit > > const &hits)
:
double Vy(const int i=0) const
Encapsulate the geometry of one entire detector (near, far, ndos)
Atom< double > MinPhysicsFrac
art::Ptr< simb::MCTruth > neutrinoInt
Truth information about neutrino interaction.
std::vector< const sim::Particle * > HitsToParticle(const std::vector< const rb::CellHit * > &hits) const
Returns vector of sim::Particle objects contributing to the given collection of hits.
std::vector< const sim::Particle * > MCTruthToParticles(art::Ptr< simb::MCTruth > const &mct) const
double purity
Purity (based on FLS energy) of neutrino interaction relative to slice.
int EveId(const int trackID) const
bool CloseInTime(const rb::CellHit &chit, const sim::PhotonSignal &phot) const
Internal helper, decides if a PhotonSignal and CellHit are close enough to match. ...