42 class ClusterTrackAna;
74 std::array<float, 5>
Cnts{{0}};
77 std::array<float, 5>
ESums{{0}};
79 std::array<float, 5>
PSums{{0}};
92 unsigned int& firstHitIndex,
93 unsigned int& lastHitIndex);
101 bool validModuleLabel =
false;
109 <<
"You must specify either a ClusterModuleLabel OR a TrackModuleLabel\n";
112 int tmp = pset.get<
int>(
"TruthOrigin", 0);
115 if (pset.has_key(
"SkipPDGCodes"))
fSkipPDGCodes = pset.get<std::vector<int>>(
"SkipPDGCodes");
116 fBadEP = pset.get<
float>(
"BadEP", 0.);
118 int intpc = pset.get<
int>(
"InTPC", -1);
119 if (intpc >= 0)
fInTPC = intpc;
134 auto const* geom = lar::providerFrom<geo::Geometry>();
143 <<
"Failed to get a handle to MCParticles using largeant\n";
161 std::vector<int> trkIDs;
163 std::vector<unsigned int> MCPIs;
170 bool isCharged = (pdg == 11) || (pdg == 13) || (pdg == 211) || (pdg == 321) || (pdg == 2212);
171 if (!isCharged)
continue;
173 float TMeV = 1000 * (part->
E() - part->
Mass());
175 if (TMeV < 1)
continue;
177 unsigned int mcpi = UINT_MAX;
178 for (mcpi = 0; mcpi < (*mcpHandle).size(); ++mcpi)
179 if ((*mcpHandle)[mcpi].TrackId() == trackID)
break;
180 if (mcpi == UINT_MAX) {
181 std::cout <<
"Failed to find a MCParticle from ParticleList\n";
184 trkIDs.push_back(trackID);
185 MCPIs.push_back(mcpi);
187 if (trkIDs.empty())
return;
193 std::vector<std::pair<unsigned int, unsigned int>> hitRange(geom->NTPC() + 1);
197 for (
auto& hr : hitRange)
198 hr = std::make_pair(UINT_MAX, UINT_MAX);
199 for (
size_t iht = 0; iht < (*fHitHandle).size(); ++iht) {
200 auto&
hit = (*fHitHandle)[iht];
202 unsigned int tpc =
hit.WireID().TPC;
210 for (
auto& tide : tides) {
212 if (tide.energyFrac < 0.5)
continue;
213 int trackID = tide.trackID;
215 for (
size_t indx = 0; indx < trkIDs.size(); ++indx) {
216 if (trkIDs[indx] == trackID) {
226 if (hitRange[tpc].first == UINT_MAX) hitRange[tpc].first = iht;
227 hitRange[tpc].second = iht;
238 myprt <<
"Checking " << trkIDs.size() <<
" selected MCParticles with the requested TruthOrigin";
239 if (
fInTPC == UINT_MAX) { myprt <<
" in all TPCs\n"; }
241 myprt <<
" in TPC " <<
fInTPC;
243 myprt <<
" in Run " << evt.
run() <<
" Event " << evt.
event();
245 myprt <<
"Found " << nMatch <<
" MC-matched hits with the requested TruthOrigin";
246 myprt <<
" out of " << nInTPC <<
" total hits.\n";
247 myprt <<
"Found " << noMatch <<
" hits not matched to ANY MCParticle.\n";
250 std::vector<std::pair<unsigned int, unsigned int>> mcpCnt;
252 if (mcpi == UINT_MAX)
continue;
253 unsigned short mIndx = 0;
254 for (mIndx = 0; mIndx < mcpCnt.size(); ++mIndx)
255 if (mcpCnt[mIndx].first == mcpi)
break;
256 if (mIndx == mcpCnt.size()) mcpCnt.push_back(std::make_pair(mcpi, 0));
258 ++mcpCnt[mIndx].second;
260 myprt <<
" MCPI TrackID PDGCode T(MeV) nHits Process";
261 for (
auto mcpcnt : mcpCnt) {
263 unsigned int mcpi = mcpcnt.first;
264 auto& mcp = (*mcpHandle)[mcpi];
265 myprt << std::setw(5) << mcpi;
266 myprt << std::setw(8) << mcp.TrackId();
267 myprt << std::setw(8) <<
abs(mcp.PdgCode());
268 int TMeV = 1000 * (mcp.E() - mcp.Mass());
269 myprt << std::setw(9) << TMeV;
270 myprt << std::setw(8) << mcpcnt.second;
271 myprt <<
" " << mcp.Process();
276 std::vector<std::vector<unsigned int>> recoHits;
278 std::vector<unsigned int> recoIndex;
290 if (!hitsFromCls.isValid())
291 throw cet::exception(
"ClusterTrackAna") <<
"Failed to get a handle to Cluster -> Hit assns\n";
292 for (
unsigned int icl = 0; icl < (*inputClusters).size(); ++icl) {
293 std::vector<art::Ptr<recob::Hit>> cluhits = hitsFromCls.at(icl);
294 if (cluhits.empty())
continue;
298 <<
"Hits associated with ClusterModuleLabel are in a different collection than " 305 std::vector<unsigned int> hitsIndex;
306 for (
auto& cluhit : cluhits) {
307 if (
fHitMCPIndex[cluhit.key()] != UINT_MAX) hitsIndex.push_back(cluhit.key());
309 if (hitsIndex.empty())
continue;
310 recoIndex.push_back(icl);
311 recoHits.push_back(hitsIndex);
320 if (!hitsFromTrk.isValid())
321 throw cet::exception(
"ClusterTrackAna") <<
"Failed to get a handle to Track -> Hit assns\n";
322 for (
unsigned int itk = 0; itk < (*inputTracks).size(); ++itk) {
323 std::vector<art::Ptr<recob::Hit>> trkhits = hitsFromTrk.at(itk);
324 if (trkhits.empty())
continue;
328 <<
"Hits associated with TrackModuleLabel are in a different collection than " 332 std::vector<unsigned int> hitsIndex;
333 for (
auto& trkhit : trkhits) {
334 if (
fHitMCPIndex[trkhit.key()] != UINT_MAX) hitsIndex.push_back(trkhit.key());
336 if (hitsIndex.empty())
continue;
337 recoIndex.push_back(itk);
338 recoHits.push_back(hitsIndex);
342 for (
const auto& tpcid : geom->Iterate<
geo::TPCID>()) {
343 unsigned int tpc = tpcid.TPC;
344 if (hitRange[tpc].first == UINT_MAX)
continue;
346 for (
unsigned short plane = 0; plane < geom->Nplanes(); ++plane) {
347 unsigned int tpcMatHitCnt = 0;
348 unsigned int tpcTotHitCnt = 0;
351 std::vector<std::pair<unsigned int, float>> mcpCnt;
353 std::vector<std::vector<std::pair<unsigned int, float>>> mcpClsCnt;
354 for (
unsigned int iht = hitRange[tpc].first; iht <= hitRange[tpc].second; ++iht) {
355 auto&
hit = (*fHitHandle)[iht];
356 if (
hit.WireID().TPC != tpc)
continue;
357 if (
hit.WireID().Plane != plane)
continue;
364 unsigned short mIndx = 0;
365 for (mIndx = 0; mIndx < mcpCnt.size(); ++mIndx)
366 if (mcpCnt[mIndx].first == mcpi)
break;
367 if (mIndx == mcpCnt.size()) {
369 mcpCnt.push_back(std::make_pair(mcpi, 0));
370 mcpClsCnt.resize(mcpCnt.size());
373 ++mcpCnt[mIndx].second;
378 myprt <<
"TPC:" << tpc <<
" Plane:" << plane <<
" has " << tpcMatHitCnt <<
"/" 380 myprt <<
" (MC-matched hits) / (all hits)";
382 if (tpcMatHitCnt < 3)
continue;
384 for (
unsigned int ii = 0; ii < recoHits.size(); ++ii) {
385 float nRecoHitsInPlane = 0;
386 float nRecoMatHitsInPlane = 0;
387 for (
auto iht : recoHits[ii]) {
388 auto&
hit = (*fHitHandle)[iht];
389 if (
hit.WireID().TPC != tpc)
continue;
390 if (
hit.WireID().Plane != plane)
continue;
393 ++nRecoMatHitsInPlane;
397 unsigned short mIndx = 0;
398 for (mIndx = 0; mIndx < mcpCnt.size(); ++mIndx)
399 if (mcpCnt[mIndx].first == mcpi)
break;
400 if (mIndx == mcpCnt.size()) {
401 std::cout <<
"Logic error: fHitMCPIndex = " <<
fHitMCPIndex[iht]
402 <<
" is valid but isn't in the list of MCParticles to use. Please send email " 403 "to baller@fnal.gov.\n";
406 unsigned short cIndx = 0;
407 for (cIndx = 0; cIndx < mcpClsCnt[mIndx].size(); ++cIndx)
408 if (mcpClsCnt[mIndx][cIndx].first == ii)
break;
409 if (cIndx == mcpClsCnt[mIndx].
size()) mcpClsCnt[mIndx].push_back(std::make_pair(ii, 0));
410 ++mcpClsCnt[mIndx][cIndx].second;
412 if (nRecoMatHitsInPlane == 0)
continue;
415 for (
unsigned short mIndx = 0; mIndx < mcpCnt.size(); ++mIndx) {
417 if (mcpCnt[mIndx].
second < 3)
continue;
418 unsigned int mcpi = mcpCnt[mIndx].first;
419 auto& mcp = (*mcpHandle)[mcpi];
420 unsigned int pdgCode =
abs(mcp.PdgCode());
421 unsigned short pIndx = USHRT_MAX;
422 if (pdgCode == 11) pIndx = 0;
423 if (pdgCode == 13) pIndx = 1;
424 if (pdgCode == 211) pIndx = 2;
425 if (pdgCode == 321) pIndx = 3;
426 if (pdgCode == 2212) pIndx = 4;
427 if (mcpClsCnt[mIndx].
empty()) {
433 myprt <<
" MCPI " << mcpi;
434 int TMeV = 1000 * (mcp.E() - mcp.Mass());
435 myprt <<
" " << TMeV <<
" MeV";
437 if (pdgCode == 11) partName =
"El";
438 if (pdgCode == 13) partName =
"Mu";
439 if (pdgCode == 211) partName =
"Pi";
440 if (pdgCode == 311) partName =
"Ka";
441 if (pdgCode == 2212) partName =
"Pr";
442 myprt << std::setw(5) << partName;
444 unsigned int firstHitIndex = UINT_MAX;
445 unsigned int lastHitIndex = UINT_MAX;
447 myprt <<
" Failed to reconstruct. Truth-matched hit range from ";
451 myprt <<
" <- EP = 0!";
455 std::pair<unsigned int, float> big = std::make_pair(UINT_MAX, 0);
456 for (
unsigned short cIndx = 0; cIndx < mcpClsCnt[mIndx].size(); ++cIndx) {
457 auto& mcc = mcpClsCnt[mIndx][cIndx];
458 if (mcc.second > big.second) big = mcc;
460 if (big.first == UINT_MAX) {
463 unsigned int mcpi = mcpCnt[mIndx].first;
464 auto& mcp = (*mcpHandle)[mcpi];
465 myprt <<
"Match failed: mcpi " << mcpi <<
" pdg " << mcp.PdgCode();
467 std::cout <<
"match failed mIndx " << mIndx <<
"\n";
470 unsigned int ii = big.first;
471 float eff = big.second / mcpCnt[mIndx].second;
472 float nRecoHitsInPlane = 0;
474 unsigned int firstRecoHitIndex = UINT_MAX;
475 unsigned int lastRecoHitIndex = UINT_MAX;
476 for (
auto iht : recoHits[ii]) {
477 auto&
hit = (*fHitHandle)[iht];
478 if (
hit.WireID().TPC != tpc)
continue;
479 if (
hit.WireID().Plane != plane)
continue;
481 if (firstRecoHitIndex == UINT_MAX) {
482 firstRecoHitIndex = iht;
483 lastRecoHitIndex = iht;
485 unsigned int wire = (*fHitHandle)[iht].WireID().Wire;
486 if (wire < (*
fHitHandle)[firstRecoHitIndex].
WireID().Wire) firstRecoHitIndex = iht;
487 if (wire > (*
fHitHandle)[lastRecoHitIndex].
WireID().Wire) lastRecoHitIndex = iht;
489 float pur = big.second / nRecoHitsInPlane;
491 float ep = eff * pur;
495 bool hasBadEP = (ep <
fBadEP);
501 myprt <<
"Hit location format is TPC:Plane:Wire:Tick\n";
502 myprt <<
" MCPI Ptcl T(MeV) nMCPHits ____From____ _____To_____";
503 if (inputClusters.
isValid()) { myprt <<
" clsID"; }
508 <<
" ____From____ _____To_____ nRecoHits nMCPRecoHits Eff Pur EP\n";
511 myprt << std::setw(5) << mcpi;
514 if (pdgCode == 11) partName =
" Electron";
515 if (pdgCode == 13) partName =
" Muon";
516 if (pdgCode == 211) partName =
" Pion";
517 if (pdgCode == 311) partName =
" Kaon";
518 if (pdgCode == 2212) partName =
" Proton";
520 int TMeV = 1000 * (mcp.E() - mcp.Mass());
521 myprt << std::setw(7) << TMeV;
523 myprt << std::setw(10) << mcpCnt[mIndx].second;
524 unsigned int firstTruHitIndex = UINT_MAX;
525 unsigned int lastTruHitIndex = UINT_MAX;
527 myprt << std::setw(14) <<
HitLocation(firstTruHitIndex);
528 myprt << std::setw(14) <<
HitLocation(lastTruHitIndex);
532 auto& cls = (*inputClusters)[recoIndex[ii]];
535 else if (inputTracks.
isValid()) {
537 auto&
trk = (*inputTracks)[recoIndex[ii]];
540 myprt << std::setw(6) << id;
541 myprt << std::setw(14) <<
HitLocation(firstRecoHitIndex);
542 myprt << std::setw(14) <<
HitLocation(lastRecoHitIndex);
543 myprt << std::setw(12) << (int)nRecoHitsInPlane;
544 myprt << std::setw(13) << (int)big.second;
545 myprt << std::fixed << std::setprecision(2);
546 myprt << std::setw(8) << eff;
547 myprt << std::setw(8) << pur;
548 myprt << std::setw(8) << ep;
549 if (hasBadEP) myprt <<
" <- BadEP";
550 myprt <<
" Evt: " << evt.
event();
565 if (iht >= (*fHitHandle).size())
return "NA";
566 auto&
hit = (*fHitHandle)[iht];
575 unsigned int& firstHitIndex,
576 unsigned int& lastHitIndex)
580 firstHitIndex = UINT_MAX;
581 lastHitIndex = UINT_MAX;
582 for (
unsigned int iht = 0; iht < (*fHitHandle).size(); ++iht) {
584 auto&
hit = (*fHitHandle)[iht];
585 if (
hit.WireID().TPC != tpc)
continue;
586 if (
hit.WireID().Plane != plane)
continue;
587 if (firstHitIndex == UINT_MAX) {
591 unsigned int wire = (*fHitHandle)[iht].WireID().Wire;
602 myprt <<
"ClusterTrackAna summary results for " <<
fEventCnt;
610 if (
fInTPC != UINT_MAX) { myprt <<
" in TPC " <<
fInTPC; }
612 myprt <<
" in all TPCs";
616 for (
unsigned short pIndx = 0; pIndx < 5; ++pIndx)
619 myprt <<
"No ClusterTrackAna results";
625 myprt <<
"Efficiency (Eff), Purity (Pur) and Eff * Pur (EP) by selected truth particle types\n";
626 std::array<std::string, 5> pName = {{
"El",
"Mu",
"Pi",
"K ",
"P "}};
627 myprt <<
"particle Eff Pur EP\n";
628 myprt <<
"-------------------------------";
629 for (
unsigned short pIndx = 0; pIndx < 5; ++pIndx) {
630 if (
Cnts[pIndx] == 0)
continue;
632 myprt <<
"\n " << pName[pIndx] <<
" ";
633 myprt << std::fixed << std::setprecision(3);
635 myprt << std::setw(8) << ave;
637 myprt << std::setw(8) << ave;
639 myprt << std::setw(8) << ave;
640 if (pIndx == 0)
continue;
642 sumE +=
ESums[pIndx];
643 sumP +=
PSums[pIndx];
645 if (cnts == 0)
return;
647 myprt <<
"Averages for all selected truth particles\n";
648 myprt <<
"Ave Eff " << sumE / cnts;
649 myprt <<
" Ave Pur " << sumP / cnts;
650 myprt <<
" Ave EP " << sumEP / cnts;
651 myprt <<
" nBadEP " <<
fNBadEP;
652 myprt <<
" (EP < " << std::fixed << std::setprecision(2) <<
fBadEP <<
")";
654 myprt <<
"MCParticle counts in all TPCs and Planes:";
655 for (
unsigned short pIndx = 0; pIndx < 5; ++pIndx) {
656 if (
Cnts[pIndx] == 0)
continue;
657 myprt <<
" " << pName[pIndx] <<
" " << (int)
Cnts[pIndx];
double E(const int i=0) const
MaybeLogger_< ELseverityLevel::ELsev_info, true > LogVerbatim
art::InputTag fTrackModuleLabel
ClusterTrackAna(fhicl::ParameterSet const &p)
Utilities related to art service access.
std::vector< sim::TrackIDE > HitToTrackIDEs(detinfo::DetectorClocksData const &clockData, recob::Hit const &hit) const
art::InputTag fClusterModuleLabel
ClusterTrackAna & operator=(ClusterTrackAna const &)=delete
Declaration of signal hit object.
std::array< float, 5 > PSums
list_type::const_iterator const_iterator
simb::Origin_t Origin() const
enum simb::_ev_origin Origin_t
event origin types
constexpr auto abs(T v)
Returns the absolute value of the argument.
Cluster finding and building.
std::string HitLocation(unsigned int iht)
packs hit WireID and PeakTime into a compact format
std::array< float, 5 > Cnts
EDAnalyzer(fhicl::ParameterSet const &pset)
bool isValid() const noexcept
decltype(auto) constexpr size(T &&obj)
ADL-aware version of std::size.
void analyze(art::Event const &e) override
#define DEFINE_ART_MODULE(klass)
IDparameter< geo::WireID > WireID
Member type of validated geo::WireID parameter.
decltype(auto) constexpr to_string(T &&obj)
ADL-aware version of std::to_string.
const art::Ptr< simb::MCTruth > & TrackIdToMCTruth_P(int id) const
std::array< float, 5 > EPSums
Provides recob::Track data product.
EventNumber_t event() const
void FirstLastHitInPlane(unsigned int tpc, unsigned int plane, unsigned int mcpi, unsigned int &firstHitIndex, unsigned int &lastHitIndex)
The data type to uniquely identify a TPC.
Declaration of cluster object.
Detector simulation of raw signals on wires.
const sim::ParticleList & ParticleList() const
std::vector< int > fSkipPDGCodes
bool getByLabel(std::string const &label, std::string const &instance, Handle< PROD > &result) const
art::InputTag fHitModuleLabel
simb::Origin_t fTruthOrigin
bool fCompareProductIDs
compare Hit and Cluster-> Hit art product IDs on the first event
art::Handle< std::vector< recob::Hit > > fHitHandle
std::array< float, 5 > ESums
second_as<> second
Type of time stored in seconds, in double precision.
art framework interface to geometry description
cet::coded_exception< error, detail::translate > exception
decltype(auto) constexpr empty(T &&obj)
ADL-aware version of std::empty.
std::vector< unsigned int > fHitMCPIndex