diff --git a/fcl/prolog.fcl b/fcl/prolog.fcl index 9d072e8..cd1ca37 100644 --- a/fcl/prolog.fcl +++ b/fcl/prolog.fcl @@ -97,13 +97,15 @@ dioLLWeight : { } genCountLogger: { module_type: GenEventCountReader } -AllOpt : { fillMC : true +AllOpt : { + fillMC : true trkqual : "TrkQual" fillTrkQual : true trkpid : "TrkPID" fillTrkPID : true fillHits : true genealogyDepth : -1 + matchDepth : -1 } DeM : { input : "KK" diff --git a/fcl/prolog_trigger.fcl b/fcl/prolog_trigger.fcl index cf1ae1d..721e1e1 100644 --- a/fcl/prolog_trigger.fcl +++ b/fcl/prolog_trigger.fcl @@ -6,13 +6,6 @@ BEGIN_PROLOG TrkAnaTrigger : { TrkAnaTT : { @table::TrkAnaTreeMaker - candidate : { - @table::DeM - options : { - @table::AllOpt - fillHits : true - } - } supplements : [] ExtraMCStepCollectionTags : [] diagLevel : 2 @@ -37,5 +30,33 @@ TrkAnaTrigger : { } } } +TrkAnaTrigger : { + @table::TrkAnaTrigger + TrkAnaDeMTT : { + @table::TrkAnaTrigger.TrkAnaTT + candidate : { + branch: "trk" + suffix : "DeM" + options : { + @table::AllOpt + fillHits : true + } + } + } + + TrkAnaDePTT : { + @table::TrkAnaTrigger.TrkAnaTT + candidate : { + branch: "trk" + suffix : "DeP" + options : { + @table::AllOpt + fillHits : true + } + } + } + + +} END_PROLOG diff --git a/inc/InfoMCStructHelper.hh b/inc/InfoMCStructHelper.hh index e0afdc0..ac0f608 100644 --- a/inc/InfoMCStructHelper.hh +++ b/inc/InfoMCStructHelper.hh @@ -60,7 +60,7 @@ namespace mu2e { void fillTrkInfoMC(const KalSeed& kseed, const KalSeedMC& kseedmc, TrkInfoMC& trkinfomc); void fillTrkInfoMCDigis(const KalSeed& kseed, const KalSeedMC& kseedmc, TrkInfoMC& trkinfomc); void fillHitInfoMC(const KalSeedMC& kseedmc, TrkStrawHitInfoMC& tshinfomc, const TrkStrawHitMC& tshmc); - void fillAllSimInfos(const KalSeedMC& kseedmc, const PrimaryParticle& primary, std::vector& siminfos, int n_generations); + void fillAllSimInfos(const KalSeedMC& kseedmc, const PrimaryParticle& primary, std::vector& siminfos, int n_generations, int n_match); void fillVDInfo(KalSeed const& kseed, const KalSeedMC& kseedmc, std::vector& vdinfos); void fillHitInfoMCs(const KalSeedMC& kseedmc, std::vector& tshinfomcs); void fillCaloClusterInfoMC(CaloClusterMC const& ccmc, CaloClusterInfoMC& ccimc); diff --git a/inc/SimInfo.hh b/inc/SimInfo.hh index 2ded05a..0615e42 100644 --- a/inc/SimInfo.hh +++ b/inc/SimInfo.hh @@ -13,13 +13,18 @@ namespace mu2e // general info about the SimParticle which was simulated struct SimInfo { bool valid = false; // true/false as to whether the data is valid + Int_t nhits = -1; // number of matched hits in the track generated by this particle + Int_t nactive = -1; // number of active matched hits in the track generated by this particle + Int_t rank = -1; // rank of this particle, by # of associated hits Int_t pdg = -1; // true PDG code + Int_t proc = -1; // process code Int_t gen = -1; // generator code - Float_t time = -1.0; // time of this step - XYZVectorF mom = XYZVectorF(); // particle momentum at the start of this step - XYZVectorF pos = XYZVectorF(); // particle position at the start of this step + Float_t time = -1.0; // Origin time of the SimParticle + Int_t index; // index into the SimInfo vector + XYZVectorF mom = XYZVectorF(); // origin momentumof the SimParticle + XYZVectorF pos = XYZVectorF(); // origin position of the SimParticle MCRelationship prirel = MCRelationship::none; // relationship to the event primary particles - MCRelationship trkrel = MCRelationship::none; // relationship to the particle that created the track + MCRelationship trkrel = MCRelationship::none; // relationship to the particle that created hits in the track void reset() { *this = SimInfo(); } }; } diff --git a/src/InfoMCStructHelper.cc b/src/InfoMCStructHelper.cc index 3f29ac7..aee599d 100644 --- a/src/InfoMCStructHelper.cc +++ b/src/InfoMCStructHelper.cc @@ -123,39 +123,51 @@ namespace mu2e { tshinfomc.doca = -1*dperp; } - void InfoMCStructHelper::fillAllSimInfos(const KalSeedMC& kseedmc, const PrimaryParticle& primary, std::vector& siminfos, int n_generations) { - auto trkprimaryptr = kseedmc.simParticle().simParticle(_spcH); - auto trkprimary = trkprimaryptr->originParticle(); - - auto current_sim_particle_ptr = trkprimaryptr; - auto current_sim_particle = trkprimary; - if (n_generations == -1) { // means do all generations + void InfoMCStructHelper::fillAllSimInfos(const KalSeedMC& kseedmc, const PrimaryParticle& primary, std::vector& siminfos, int n_generations, int n_match) { + // interpret -1 as no llimit + if (n_generations == -1) { n_generations = std::numeric_limits::max(); } - - for (int i_generation = 0; i_generation < n_generations; ++i_generation) { - SimInfo sim_info; - fillSimInfo(current_sim_particle, sim_info); - sim_info.trkrel = MCRelationship(current_sim_particle_ptr, trkprimaryptr); - - auto bestprimarysp = primary.primarySimParticles().front(); - MCRelationship bestrel; - for(auto const& spp : primary.primarySimParticles()){ - MCRelationship mcrel(current_sim_particle_ptr, spp); - if(mcrel > bestrel){ - bestrel = mcrel; - bestprimarysp = spp; + if (n_match == -1 ) { + n_match = std::numeric_limits::max(); + } + for(int imatch = 0 ; imatch < std::min(n_match,static_cast(kseedmc.simParticles().size())); ++imatch) { + auto trkprimaryptr = kseedmc.simParticle(imatch).simParticle(_spcH); + auto trkprimary = trkprimaryptr->originParticle(); + auto current_sim_particle_ptr = trkprimaryptr; + auto current_sim_particle = trkprimary; + + for (int i_generation = 0; i_generation < n_generations; ++i_generation) { + SimInfo sim_info; + fillSimInfo(current_sim_particle, sim_info); + sim_info.trkrel = MCRelationship(current_sim_particle_ptr, trkprimaryptr); + sim_info.rank = imatch; + + auto bestprimarysp = primary.primarySimParticles().front(); + MCRelationship bestrel; + for(auto const& spp : primary.primarySimParticles()){ + MCRelationship mcrel(current_sim_particle_ptr, spp); + if(mcrel > bestrel){ + bestrel = mcrel; + bestprimarysp = spp; + } + } + sim_info.prirel = bestrel; + // only count hits for direct contributors + if(i_generation == 0){ + sim_info.nhits = kseedmc.simParticle(imatch)._nhits; + sim_info.nactive = kseedmc.simParticle(imatch)._nactive; + } + // record the index this object will have + sim_info.index = siminfos.size(); + siminfos.push_back(sim_info); + if (current_sim_particle.parent().isNonnull()) { + current_sim_particle_ptr = current_sim_particle.parent(); + current_sim_particle = current_sim_particle_ptr->originParticle(); + } + else { + break; // this particle doesn't have a parent } - } - sim_info.prirel = bestrel; - - siminfos.push_back(sim_info); - if (current_sim_particle.parent().isNonnull()) { - current_sim_particle_ptr = current_sim_particle.parent(); - current_sim_particle = current_sim_particle_ptr->originParticle(); - } - else { - break; // this particle doesn't have a parent } } @@ -173,12 +185,13 @@ namespace mu2e { } } if (!already_added) { + auto trkprimaryptr = kseedmc.simParticle().simParticle(_spcH); sim_info.trkrel = MCRelationship(spp, trkprimaryptr); sim_info.prirel = MCRelationship(spp, spp); + sim_info.index = siminfos.size(); siminfos.push_back(sim_info); } } - } @@ -192,12 +205,11 @@ namespace mu2e { } void InfoMCStructHelper::fillSimInfo(const SimParticle& sp, SimInfo& siminfo) { - GeomHandle det; - siminfo.valid = true; + if(sp.genParticle().isNonnull())siminfo.gen = sp.genParticle()->generatorId().id(); + siminfo.proc = sp.creationCode(); siminfo.pdg = sp.pdgId(); - siminfo.gen = sp.creationCode(); siminfo.time = sp.startGlobalTime(); siminfo.mom = XYZVectorF(sp.startMomentum()); siminfo.pos = XYZVectorF(det->toDetector(sp.startPosition())); diff --git a/src/TrkAnaTreeMaker_module.cc b/src/TrkAnaTreeMaker_module.cc index 58a5df3..e572864 100644 --- a/src/TrkAnaTreeMaker_module.cc +++ b/src/TrkAnaTreeMaker_module.cc @@ -113,6 +113,7 @@ namespace mu2e { fhicl::Atom filltrkpid{Name("fillTrkPID"), Comment("Switch to turn on filling of the full TrkPIDInfo for this set of tracks"), false}; fhicl::Atom required{Name("required"), Comment("True/false if you require this type of track in the event"), false}; fhicl::Atom genealogyDepth{Name("genealogyDepth"), Comment("The depth of the genealogy information you want to keep"), 1}; + fhicl::Atom matchDepth{Name("matchDepth"), Comment("The depth into the MC true particle matching you want to keep"), 1}; }; struct BranchConfig { @@ -854,7 +855,7 @@ namespace mu2e { _infoMCStructHelper.fillTrkInfoMC(kseed, kseedmc, _allMCTIs.at(i_branch)); auto& mcvdis = _allMCVDInfos.at(i_branch); _infoMCStructHelper.fillVDInfo(kseed, kseedmc, mcvdis); - _infoMCStructHelper.fillAllSimInfos(kseedmc, primary, _allMCSimTIs.at(i_branch), branchConfig.options().genealogyDepth()); + _infoMCStructHelper.fillAllSimInfos(kseedmc, primary, _allMCSimTIs.at(i_branch), branchConfig.options().genealogyDepth(), branchConfig.options().matchDepth()); if(_conf.diag() > 1 || (_conf.fillhits() && branchConfig.options().fillhits())){ _infoMCStructHelper.fillHitInfoMCs(kseedmc, _allTSHIMCs.at(i_branch));