From 5e07bde056fc27e5d553efb22754dd58d3df37d2 Mon Sep 17 00:00:00 2001 From: Andrew Edmonds Date: Fri, 16 Dec 2022 12:53:13 -0600 Subject: [PATCH] Add ability to configure virtual detectors in TrkAna and include an example --- fcl/TrkAnaReco_addStopTgtVDs.fcl | 19 +++++++ fcl/prolog.fcl | 14 +++++ inc/InfoMCStructHelper.hh | 2 +- src/InfoMCStructHelper.cc | 2 +- src/TrkAnaTreeMaker_module.cc | 89 ++++++++++++++++++-------------- 5 files changed, 86 insertions(+), 40 deletions(-) create mode 100644 fcl/TrkAnaReco_addStopTgtVDs.fcl diff --git a/fcl/TrkAnaReco_addStopTgtVDs.fcl b/fcl/TrkAnaReco_addStopTgtVDs.fcl new file mode 100644 index 0000000..3cea34d --- /dev/null +++ b/fcl/TrkAnaReco_addStopTgtVDs.fcl @@ -0,0 +1,19 @@ +// This fcl file gives an example of how an analyzer +// can write a TrkAna tree with different VirtualDetectors +// +// This is useful if .... +// +// In this example, we will write branches for ST_Out virtual detector +// + +// This example builds on TrkAnaReco.fcl +#include "TrkAna/fcl/TrkAnaReco.fcl" + +// First define the suffix of the branches we want +// In this example, "tgt" will mean we create branched "detgt", "demctgt", "uetgt", etc. +physics.analyzers.TrkAnaNeg.candidate.segmentSuffixes : [ @sequence::physics.analyzers.TrkAnaNeg.candidate.segmentSuffixes, "tgt" ] + +// Then define the virtual detectors we want to be included in this branch +// NB we add a sequence because we sometimes combine virtual detectors for one "interesting" plane +// (e.g. TT_FrontHollow and TT_FrontPA define the entrance to the tracker) +physics.analyzers.TrkAnaNeg.candidate.segmentVIDs : [ @sequence::physics.analyzers.TrkAnaNeg.candidate.segmentVIDs, ["ST_Out"] ] diff --git a/fcl/prolog.fcl b/fcl/prolog.fcl index 6720fe8..5d52fbc 100644 --- a/fcl/prolog.fcl +++ b/fcl/prolog.fcl @@ -133,6 +133,11 @@ dioLLWeight : { genCountLogger: { module_type: GenEventCountReader } # candidate configuration for TrkAna +StdSegments : { + segmentSuffixes : [ "ent", "mid", "xit" ] + segmentVIDs : [ [ "TT_FrontHollow", "TT_FrontPA"], [ "TT_Mid", "TT_MidInner"], ["TT_Back"] ] +} + AllOpt : { fillMC : true trkqual : "TrkQual" fillTrkQual : true @@ -141,6 +146,7 @@ AllOpt : { fillMC : true fillHits : true genealogyDepth : 5 fillBestCrv : true + bestCrvModules : [ "BestCrv" ] bestCrvInstances : [ "first" ] bestCrvBranches : [ "bestcrv" ] @@ -149,34 +155,42 @@ AllOpt : { fillMC : true DeM : { input : "KFF" branch : "de" suffix : "DeM" + @table::StdSegments } UeM : { input : "KFF" branch : "ue" suffix : "UeM" + @table::StdSegments } DmuM : { input : "KFF" branch : "dm" suffix : "DmuM" + @table::StdSegments } UmuM : { input : "KFF" branch : "um" suffix : "UmuM" + @table::StdSegments } DeP : { input : "KFF" branch : "de" suffix : "DeP" + @table::StdSegments } UeP : { input : "KFF" branch : "ue" suffix : "UeP" + @table::StdSegments } DmuP : { input : "KFF" branch : "dm" suffix : "DmuP" + @table::StdSegments } UmuP : { input : "KFF" branch : "um" suffix : "UmuP" + @table::StdSegments } diff --git a/inc/InfoMCStructHelper.hh b/inc/InfoMCStructHelper.hh index f5014b8..b72e6d3 100644 --- a/inc/InfoMCStructHelper.hh +++ b/inc/InfoMCStructHelper.hh @@ -57,7 +57,7 @@ namespace mu2e { void fillHitInfoMC(const KalSeedMC& kseedmc, TrkStrawHitInfoMC& tshinfomc, const TrkStrawHitMC& tshmc); void fillAllSimInfos(const KalSeedMC& kseedmc, std::vector& siminfos, int n_generations); void fillPriInfo(const KalSeedMC& kseedmc, const PrimaryParticle& primary, SimInfo& priinfo); - void fillTrkInfoMCStep(const KalSeedMC& kseedmc, TrkInfoMCStep& trkinfomcstep, std::vector const& vids, double target_time); + void fillTrkInfoMCStep(const KalSeedMC& kseedmc, TrkInfoMCStep& trkinfomcstep, std::vector const& vids, double target_time); void fillHitInfoMCs(const KalSeedMC& kseedmc, std::vector& tshinfomcs); void fillCaloClusterInfoMC(CaloClusterMC const& ccmc, CaloClusterInfoMC& ccimc); diff --git a/src/InfoMCStructHelper.cc b/src/InfoMCStructHelper.cc index 946dc98..6dbbc85 100644 --- a/src/InfoMCStructHelper.cc +++ b/src/InfoMCStructHelper.cc @@ -164,7 +164,7 @@ namespace mu2e { } void InfoMCStructHelper::fillTrkInfoMCStep(const KalSeedMC& kseedmc, TrkInfoMCStep& trkinfomcstep, - std::vector const& vids, double target_time) { + std::vector const& vids, double target_time) { GeomHandle det; GlobalConstantsHandle pdt; diff --git a/src/TrkAnaTreeMaker_module.cc b/src/TrkAnaTreeMaker_module.cc index 58f8937..feba22c 100644 --- a/src/TrkAnaTreeMaker_module.cc +++ b/src/TrkAnaTreeMaker_module.cc @@ -28,6 +28,7 @@ #include "Offline/GeometryService/inc/VirtualDetector.hh" #include "Offline/RecoDataProducts/inc/CrvCoincidenceCluster.hh" #include "Offline/MCDataProducts/inc/CrvCoincidenceClusterMC.hh" +#include "Offline/Mu2eUtilities/inc/fromStrings.hh" // Framework includes. #include "art/Framework/Core/EDAnalyzer.h" #include "art/Framework/Principal/Event.h" @@ -123,6 +124,8 @@ namespace mu2e { fhicl::Atom input{Name("input"), Comment("KalSeedCollection input tag (use prefix if fcl parameter suffix is defined)")}; fhicl::Atom branch{Name("branch"), Comment("Name of output branch")}; fhicl::Atom suffix{Name("suffix"), Comment("Fit suffix (e.g. DeM)"), ""}; + fhicl::Sequence segmentSuffixes{Name("segmentSuffixes"), Comment("The suffix to the branch for this segment (e.g. putting \"ent\" will give a branch \"deent\", and if fillMC = true, \"demcent\")")}; + fhicl::Sequence> segmentVIDs{Name("segmentVIDs"), Comment("The VirtualDetectorId that this segment should find tits position from")}; fhicl::Table options{Name("options"), Comment("Optional arguments for a branch")}; }; @@ -208,7 +211,9 @@ namespace mu2e { std::vector > _allKSCHs; // track branches (outputs) std::vector _allTIs; - std::vector _allEntTIs, _allMidTIs, _allXitTIs; + // std::vector _allEntTIs, _allMidTIs, _allXitTIs; + std::map> _allSegmentTIs; + std::vector _allTCHIs; // quality branches (inputs) std::vector > > _allRQCHs; // outer vector is for each candidate/supplement, inner vector is all RecoQuals @@ -229,12 +234,12 @@ namespace mu2e { art::Handle _ksmcah; art::Handle _ccmcch; art::InputTag _primaryParticleTag; - std::vector _entvids, _midvids, _xitvids; + std::map>> _allSegmentVIDs; // MC truth branches (outputs) std::vector _allMCTIs; std::vector> _allMCSimTIs; std::vector _allMCPriTIs; - std::vector _allMCEntTIs, _allMCMidTIs, _allMCXitTIs; + std::map> _allMCSegmentTIs; std::vector _allMCTCHIs; // hit level info branches @@ -315,12 +320,6 @@ namespace mu2e { _buffsize(conf().buffsize()), _splitlevel(conf().splitlevel()) { - _midvids.push_back(VirtualDetectorId::TT_Mid); - _midvids.push_back(VirtualDetectorId::TT_MidInner); - _entvids.push_back(VirtualDetectorId::TT_FrontHollow); - _entvids.push_back(VirtualDetectorId::TT_FrontPA); - _xitvids.push_back(VirtualDetectorId::TT_Back); - // collect both candidate and supplement branches into one place _allBranches.push_back(_conf.candidate()); _candidateIndex = 0; @@ -337,10 +336,24 @@ namespace mu2e { TrkInfo ti; _allTIs.push_back(ti); - TrkFitInfo ent, mid, xit; - _allEntTIs.push_back(ent); - _allMidTIs.push_back(mid); - _allXitTIs.push_back(xit); + std::vector allSegments; + std::vector allMCSegments; + const std::vector& segmentSuffixes = i_branchConfig.segmentSuffixes(); + for (size_t i_segment = 0; i_segment < segmentSuffixes.size(); ++i_segment) { + allSegments.push_back(TrkFitInfo()); // add empty structs to the vector so that ROOT can be given a location to find it + allMCSegments.push_back(TrkInfoMCStep()); + } + _allSegmentTIs[i_branch] = allSegments; + _allMCSegmentTIs[i_branch] = allMCSegments; + + const std::vector>& segmentVIDs = i_branchConfig.segmentVIDs(); + std::vector> allSegmentVIDs; + for (size_t i_segment = 0; i_segment < segmentSuffixes.size(); ++i_segment) { + std::vector thisSegmentVIDs = fromStrings(segmentVIDs.at(i_segment)); + allSegmentVIDs.push_back(thisSegmentVIDs); + } + _allSegmentVIDs[i_branch] = allSegmentVIDs; + TrkCaloHitInfo tchi; _allTCHIs.push_back(tchi); @@ -355,10 +368,6 @@ namespace mu2e { _allMCSimTIs.push_back(tempMCSimTIs); SimInfo mcpri; _allMCPriTIs.push_back(mcpri); - TrkInfoMCStep mcent, mcmid, mcxit; - _allMCEntTIs.push_back(mcent); - _allMCMidTIs.push_back(mcmid); - _allMCXitTIs.push_back(mcxit); CaloClusterInfoMC mctchi; _allMCTCHIs.push_back(mctchi); @@ -415,9 +424,12 @@ namespace mu2e { BranchConfig i_branchConfig = _allBranches.at(i_branch); std::string branch = i_branchConfig.branch(); _trkana->Branch((branch+".").c_str(),&_allTIs.at(i_branch)); - _trkana->Branch((branch+"ent.").c_str(),&_allEntTIs.at(i_branch)); - _trkana->Branch((branch+"mid.").c_str(),&_allMidTIs.at(i_branch)); - _trkana->Branch((branch+"xit.").c_str(),&_allXitTIs.at(i_branch)); + + const std::vector& segmentSuffixes = i_branchConfig.segmentSuffixes(); + for (size_t i_segment = 0; i_segment < segmentSuffixes.size(); ++i_segment) { + std::string segmentSuffix = segmentSuffixes.at(i_segment); + _trkana->Branch((branch+segmentSuffix+".").c_str(),&_allSegmentTIs.at(i_branch).at(i_segment)); + } _trkana->Branch((branch+"tch.").c_str(),&_allTCHIs.at(i_branch)); if (_conf.filltrkqual() && i_branchConfig.options().filltrkqual()) { int n_trkqual_vars = TrkQual::n_vars; @@ -482,9 +494,11 @@ namespace mu2e { _trkana->Branch((full_branchname).c_str(),&_allMCSimTIs.at(i_branch).at(i_generation)); } _trkana->Branch((branch+"mcpri.").c_str(),&_allMCPriTIs.at(i_branch),_buffsize,_splitlevel); - _trkana->Branch((branch+"mcent.").c_str(),&_allMCEntTIs.at(i_branch),_buffsize,_splitlevel); - _trkana->Branch((branch+"mcmid.").c_str(),&_allMCMidTIs.at(i_branch),_buffsize,_splitlevel); - _trkana->Branch((branch+"mcxit.").c_str(),&_allMCXitTIs.at(i_branch),_buffsize,_splitlevel); + const std::vector& segmentSuffixes = i_branchConfig.segmentSuffixes(); + for (size_t i_segment = 0; i_segment < segmentSuffixes.size(); ++i_segment) { + std::string segmentSuffix = segmentSuffixes.at(i_segment); + _trkana->Branch((branch+"mc"+segmentSuffix+".").c_str(),&_allMCSegmentTIs.at(i_branch).at(i_segment)); + } _trkana->Branch((branch+"tchmc.").c_str(),&_allMCTCHIs.at(i_branch),_buffsize,_splitlevel); // at hit-level MC information // (for the time being diagLevel will still work, but I propose removing this at some point) @@ -835,14 +849,16 @@ namespace mu2e { // get VD positions mu2e::GeomHandle vdHandle; mu2e::GeomHandle det; - const XYZVectorF& entpos = XYZVectorF(det->toDetector(vdHandle->getGlobal(*_entvids.begin()))); - const XYZVectorF& midpos = XYZVectorF(det->toDetector(vdHandle->getGlobal(*_midvids.begin()))); - const XYZVectorF& xitpos = XYZVectorF(det->toDetector(vdHandle->getGlobal(*_xitvids.begin()))); +// const XYZVectorF& entpos = XYZVectorF(det->toDetector(vdHandle->getGlobal(*_entvids.begin()))); +// const XYZVectorF& midpos = XYZVectorF(det->toDetector(vdHandle->getGlobal(*_midvids.begin()))); +// const XYZVectorF& xitpos = XYZVectorF(det->toDetector(vdHandle->getGlobal(*_xitvids.begin()))); _infoStructHelper.fillTrkInfo(kseed,_allTIs.at(i_branch)); - _infoStructHelper.fillTrkFitInfo(kseed,_allEntTIs.at(i_branch),entpos); - _infoStructHelper.fillTrkFitInfo(kseed,_allMidTIs.at(i_branch),midpos); - _infoStructHelper.fillTrkFitInfo(kseed,_allXitTIs.at(i_branch),xitpos); + const std::vector& segmentSuffixes = branchConfig.segmentSuffixes(); + for (size_t i_segment = 0; i_segment < segmentSuffixes.size(); ++i_segment) { + const XYZVectorF& pos = XYZVectorF(det->toDetector(vdHandle->getGlobal(*_allSegmentVIDs.at(i_branch).at(i_segment).begin()))); + _infoStructHelper.fillTrkFitInfo(kseed,_allSegmentTIs.at(i_branch).at(i_segment),pos); + } // _tcnt._overlaps[0] = _tcomp.nOverlap(kseed, kseed); if(_conf.diag() > 1 || (_conf.fillhits() && branchConfig.options().fillhits())){ // want hit level info @@ -925,9 +941,10 @@ namespace mu2e { auto const& kseed = *kptr; _infoMCStructHelper.fillTrkInfoMC(kseed, kseedmc, _allMCTIs.at(i_branch)); double t0 = kseed.t0().t0(); - _infoMCStructHelper.fillTrkInfoMCStep(kseedmc, _allMCEntTIs.at(i_branch), _entvids, t0); - _infoMCStructHelper.fillTrkInfoMCStep(kseedmc, _allMCMidTIs.at(i_branch), _midvids, t0); - _infoMCStructHelper.fillTrkInfoMCStep(kseedmc, _allMCXitTIs.at(i_branch), _xitvids, t0); + const std::vector& segmentSuffixes = branchConfig.segmentSuffixes(); + for (size_t i_segment = 0; i_segment < segmentSuffixes.size(); ++i_segment) { + _infoMCStructHelper.fillTrkInfoMCStep(kseedmc, _allMCSegmentTIs.at(i_branch).at(i_segment), _allSegmentVIDs.at(i_branch).at(i_segment), t0); + } _infoMCStructHelper.fillPriInfo(kseedmc, primary, _allMCPriTIs.at(i_branch)); _infoMCStructHelper.fillAllSimInfos(kseedmc, _allMCSimTIs.at(i_branch), branchConfig.options().genealogyDepth()); @@ -1006,9 +1023,7 @@ namespace mu2e { void TrkAnaTreeMaker::resetTrackBranches() { for (BranchIndex i_branch = 0; i_branch < _allBranches.size(); ++i_branch) { _allTIs.at(i_branch).reset(); - _allEntTIs.at(i_branch).reset(); - _allMidTIs.at(i_branch).reset(); - _allXitTIs.at(i_branch).reset(); + _allSegmentTIs.at(i_branch).assign(_allSegmentTIs.at(i_branch).size(), TrkFitInfo()); // we don't want to remove elements so use assign instead of clear _allTCHIs.at(i_branch).reset(); @@ -1018,9 +1033,7 @@ namespace mu2e { } _allMCPriTIs.at(i_branch).reset(); - _allMCEntTIs.at(i_branch).reset(); - _allMCMidTIs.at(i_branch).reset(); - _allMCXitTIs.at(i_branch).reset(); + _allMCSegmentTIs.at(i_branch).assign(_allMCSegmentTIs.at(i_branch).size(), TrkInfoMCStep()); _allMCTCHIs.at(i_branch).reset(); _allRQIs.at(i_branch).reset();