diff --git a/CalPatRec/fcl/prolog.fcl b/CalPatRec/fcl/prolog.fcl index 984a2d6e4f..4c21c55218 100644 --- a/CalPatRec/fcl/prolog.fcl +++ b/CalPatRec/fcl/prolog.fcl @@ -491,15 +491,18 @@ CalPatRec : { @table::CalPatRec CalTimeClusterFinder : { module_type : CalLineTimePeakFinder DiagLevel : 0 - ComboHitCollectionLabel : "makeSH" + ComboHitCollectionLabel : "makePH" CaloClusterCollectionLabel : "CaloClusterMaker" MinTimeClusterHits : 3 MinCaloClusterEnergy : 50. HitTimeSigmaThresh : 3. HitXYSigmaThresh : 3. MinHitRadius : 350. - StoppingTargetRadius : 150. + StoppingTargetRadius : 200. ClusterRadius : 50. + StrawEDepMin : 0. + StrawEDepMax : 0.003 # electron-like + IgnoreTarget : false ParticleBeta : 1. FitDirection : "downstream" } diff --git a/CalPatRec/src/CalLineTimePeakFinder_module.cc b/CalPatRec/src/CalLineTimePeakFinder_module.cc index 8322777b93..bfc9dca23f 100644 --- a/CalPatRec/src/CalLineTimePeakFinder_module.cc +++ b/CalPatRec/src/CalLineTimePeakFinder_module.cc @@ -49,6 +49,9 @@ namespace mu2e { fhicl::Atom cluster_radius {Name("ClusterRadius"), Comment("Radius around cluster position in mm for seeding")}; fhicl::Atom particle_beta {Name("ParticleBeta"), Comment("Particle beta for distance to time (0-1)"), 1.}; fhicl::Atom fit_direction {Name("FitDirection"), Comment("Fit Direction in Search (\"downstream\" or \"upstream\")") }; + fhicl::Atom hit_edep_min {Name("StrawEDepMin"), Comment("Straw hit energy minimum"), -1.}; + fhicl::Atom hit_edep_max {Name("StrawEDepMax"), Comment("Straw hit energy maximum"), -1.}; + fhicl::Atom ignore_target {Name("IgnoreTarget"), Comment("Cluster all hits in time with the calo cluster"), false}; }; //----------------------------------------------------------------------------- @@ -65,11 +68,18 @@ namespace mu2e { float cluster_radius_; // radius around cluster position in mm for seeding float particle_beta_; // particle beta (v/c) TrkFitDirection fit_dir_; // fit direction in search + float hit_edep_min_; // straw hit energy minimum + float hit_edep_max_; // straw hit energy maximum + bool ignore_target_; // ignore target origin, cluster any hits consistent in time int diag_level_; // diagnostic output //----------------------------------------------------------------------------- // Data //----------------------------------------------------------------------------- + + size_t nevents_; // counters + size_t ntime_clusters_; + const ComboHitCollection* combo_hit_col_; // input combo hit collection const CaloClusterCollection* calo_cluster_col_; // input calo cluster collection art::Handle combo_hit_col_handle_; // handle for input combo hit collection @@ -86,6 +96,7 @@ namespace mu2e { virtual void beginRun(art::Run& run ); virtual void produce (art::Event& event ); + virtual void endJob (); //----------------------------------------------------------------------------- // helper functions //----------------------------------------------------------------------------- @@ -112,7 +123,12 @@ namespace mu2e { , cluster_radius_(config().cluster_radius()) , particle_beta_(config().particle_beta()) , fit_dir_(config().fit_direction()) + , hit_edep_min_(config().hit_edep_min()) + , hit_edep_max_(config().hit_edep_max()) + , ignore_target_(config().ignore_target()) , diag_level_(config().diag_level()) + , nevents_(0) + , ntime_clusters_(0) { // declare the data products consumes(hit_tag_); @@ -152,7 +168,8 @@ namespace mu2e { evt.getByLabel(hit_tag_, combo_hit_col_handle_); if(!combo_hit_col_handle_.isValid()) { - printf("[CalLineTimePeakFinder::%s] ERROR: ComboHit collection with label \"%s\" not found! RETURN\n", __func__, hit_tag_.encode().c_str()); + printf("[%s::CalLineTimePeakFinder::%s] ERROR: ComboHit collection with label \"%s\" not found! RETURN\n", + moduleDescription().moduleLabel().c_str(), __func__, hit_tag_.encode().c_str()); combo_hit_col_ = nullptr; return false; } else { @@ -161,7 +178,8 @@ namespace mu2e { evt.getByLabel(calo_cluster_tag_, calo_cluster_col_handle_); if(!calo_cluster_col_handle_.isValid()) { - printf("[CalLineTimePeakFinder::%s] ERROR: CaloCluster collection with label \"%s\" not found! RETURN\n", __func__, calo_cluster_tag_.encode().c_str()); + printf("[%s::CalLineTimePeakFinder::%s] ERROR: CaloCluster collection with label \"%s\" not found! RETURN\n", + moduleDescription().moduleLabel().c_str(), __func__, calo_cluster_tag_.encode().c_str()); calo_cluster_col_ = nullptr; return false; } else { @@ -188,9 +206,9 @@ namespace mu2e { if(seed_dir_mag <= 0.) return; // can't define a seed direction seed_dir *= 1./seed_dir_mag; // normalize the seed direction if(std::fabs(seed_dir.z()) < 1.e-3) return; // too sharp of a slope (shouldn't be possible) - if(diag_level_ > 1) { - printf("[CalLineTimePeakFinder::%s] CaloCluster time = %.2f, position = (%.1f, %.1f, %.1f), seed direction = (%.2f, %.2f, %.2f)\n", - __func__, cl_time, cl_pos.x(), cl_pos.y(), cl_pos.z(), seed_dir.x(), seed_dir.y(), seed_dir.z()); + if(diag_level_ > 2) { + printf("[%s::CalLineTimePeakFinder::%s] CaloCluster time = %.2f, position = (%.1f, %.1f, %.1f), seed direction = (%.2f, %.2f, %.2f)\n", + moduleDescription().moduleLabel().c_str(), __func__, cl_time, cl_pos.x(), cl_pos.y(), cl_pos.z(), seed_dir.x(), seed_dir.y(), seed_dir.z()); } // Look for hits consistent this seed position and direction, and add them to the seed @@ -207,7 +225,7 @@ namespace mu2e { const double time_at_hit = cl_time + dt; // expected time at the hit z position const double time_unc = std::sqrt(hit.timeRes()*hit.timeRes() + cl.timeErr()*cl.timeErr()); // uncertainty on the time difference between the hit and the expected time based on the seed const double time_sigma = std::abs((hit_time - time_at_hit) / time_unc); // number of sigma the hit time is from the expected time based on the seed - if(diag_level_ > 2) { + if(diag_level_ > 3) { printf(" Hit %zu: time = %.2f, expected time = %.2f, time sigma = %.2f\n", i_hit, hit_time, time_at_hit, time_sigma); } @@ -216,17 +234,17 @@ namespace mu2e { // next check if the hit is consistent in space const CLHEP::Hep3Vector pos_at_hit = cl_pos + dr * seed_dir; // position along the seed direction at the same z as the hit const double cone = cluster_radius_ + std::abs(dz/dz_target) * (stopping_target_radius_ - cluster_radius_); // radius of the cone in the x-y plane at the hit z based on the seed direction and cluster position - if((pos_at_hit.perp() + cone) < min_hit_radius_) continue; // expected hit position is too low radius to still reach the tracker + if(!ignore_target_ && (pos_at_hit.perp() + cone) < min_hit_radius_) continue; // expected hit position is too low radius to still reach the tracker const double x_y_dist = (hit_pos - pos_at_hit).perp(); // distance in the x-y plane between the hit and the expected position based on the seed const double xy_unc = std::sqrt(hit.transVar() + hit.wireVar()); // uncertainty in the hit position in the x-y plane const double xy_sigma = (x_y_dist - cone) / xy_unc; // number of sigma the hit is from the cone from the cluster to the target, negative means contained within the cone - if(diag_level_ > 2) { + if(diag_level_ > 3) { printf(" Expected position at hit z: (%.1f, %.1f, %.1f), hit position: (%.1f, %.1f, %.1f), x-y distance = %.2f, target cone = %.2f, sigma = %.2f\n", pos_at_hit.x(), pos_at_hit.y(), pos_at_hit.z(), hit_pos.x(), hit_pos.y(), hit_pos.z(), x_y_dist, cone, xy_sigma); } - if(xy_sigma > hit_xy_sigma_thresh_) continue; // hit is not consistent with the seed - if(diag_level_ > 1) { - if(diag_level_ <= 2) { // print full info if not already printed + if(!ignore_target_ && xy_sigma > hit_xy_sigma_thresh_) continue; // hit is not consistent with the seed + if(diag_level_ > 2) { + if(diag_level_ <= 3) { // print full info if not already printed printf(" Hit %zu: time = %.2f, expected time = %.2f, time sigma = %.2f\n", i_hit, hit_time, time_at_hit, time_sigma); printf(" Expected position at hit z: (%.1f, %.1f, %.1f), hit position: (%.1f, %.1f, %.1f), x-y distance = %.2f, target cone = %.2f, sigma = %.2f\n", @@ -243,6 +261,7 @@ namespace mu2e { // event entry point //----------------------------------------------------------------------------- void CalLineTimePeakFinder::produce(art::Event& event ) { + ++nevents_; // output collection std::unique_ptr out_tcs(new TimeClusterCollection); @@ -261,18 +280,19 @@ namespace mu2e { const size_t n_clusters = calo_cluster_col_->size(); for(size_t i_cl = 0; i_cl < n_clusters; ++i_cl) { const auto& cl = calo_cluster_col_->at(i_cl); - if(diag_level_ > 1) { - printf("[CalLineTimePeakFinder::%s] Seeding from CaloCluster %zu with energy = %.2f, time = %.2f, N(combo hits) = %zu\n", - __func__, i_cl, cl.energyDep(), cl.time(), combo_hit_col_->size()); + if(diag_level_ > 2) { + printf("[%s::CalLineTimePeakFinder::%s] Seeding from CaloCluster %zu with energy = %.2f, time = %.2f, N(combo hits) = %zu\n", + moduleDescription().moduleLabel().c_str(), __func__, i_cl, cl.energyDep(), cl.time(), combo_hit_col_->size()); } if(cl.energyDep() < min_calo_cluster_energy_) continue; // skip clusters that don't pass the energy cut TimeCluster tc; findTimePeakInCluster(tc, cl); if(!isGoodTimeCluster(tc)) continue; // skip time clusters that don't pass selection criteria finalizeTimeCluster(tc, i_cl); - if(diag_level_ > 0) { - printf("[CalLineTimePeakFinder::%s] Found time cluster with %zu hits and t0 = %.2f\n", __func__, tc.nhits(), tc.t0().t0()); - if(diag_level_ > 1) { + if(diag_level_ > 1) { + printf("[%s::CalLineTimePeakFinder::%s] Found time cluster with %zu hits and t0 = %.1f from cluster E = %.1f MeV and t = %.1f ns\n", + moduleDescription().moduleLabel().c_str(), __func__, tc.nhits(), tc.t0().t0(), cl.energyDep(), cl.time()); + if(diag_level_ > 2) { for(size_t i = 0; i < tc.nhits(); ++i) { const size_t hit_index = tc.hits().at(i); const auto& hit = combo_hit_col_->at(hit_index); @@ -288,17 +308,26 @@ namespace mu2e { // put reconstructed time clusters into the event record //----------------------------------------------------------------------------- - if(diag_level_ > 0) { - printf("[CalLineTimePeakFinder::%s] Found %zu time clusters in total\n", __func__, out_tcs->size()); + if((diag_level_ > 0 && out_tcs->size() > 0) || (diag_level_ > 1)) { + printf("[%s::CalLineTimePeakFinder::%s] %i:%i:%i Found %zu time clusters in total\n", + moduleDescription().moduleLabel().c_str(), __func__, event.run(), event.subRun(), event.event(), out_tcs->size()); } + ntime_clusters_ += out_tcs->size(); event.put(std::move(out_tcs)); } //----------------------------------------------------------------------------- bool CalLineTimePeakFinder::isGoodHit(const ComboHit& hit) { const auto& flag = hit.flag(); - const int bkg_hit = flag.hasAnyProperty(StrawHitFlag::bkg); - if (bkg_hit != 0) return false; + // the flag being set means it was identified as background (e.g. delta electron) + const bool bkg_hit = flag.hasAnyProperty(StrawHitFlag::bkg); + // use local energy deposit testing for electron vs. proton + const float edep = hit.energyDep(); + if(diag_level_ > 3) printf(" [%s::CalLineTimePeakFinder::%s] Hit flags %10s: bkg = %o edep_flag = %o, edep = %.6f\n", + moduleDescription().moduleLabel().c_str(), __func__, flag.hex().c_str(), bkg_hit, flag.hasAnyProperty(StrawHitFlag::energysel), edep); + if(bkg_hit != 0) return false; + if(edep < hit_edep_min_) return false; + if(hit_edep_max_ > 0. && edep > hit_edep_max_) return false; return true; } @@ -333,6 +362,11 @@ namespace mu2e { tc._t0._t0err = std::sqrt(std::max(0., t0sq - avg_t0*avg_t0)); } + //----------------------------------------------------------------------------- + void CalLineTimePeakFinder::endJob() { + if(diag_level_ > 0) printf("[%s::CalLineTimePeakFinder::%s] Found %zu time clusters from %zu events\n", + moduleDescription().moduleLabel().c_str(), __func__, ntime_clusters_, nevents_); + } } // end namespace mu2e using mu2e::CalLineTimePeakFinder; diff --git a/CalPatRec/src/TZClusterFinder_module.cc b/CalPatRec/src/TZClusterFinder_module.cc index 58bbaf76bf..00d769e883 100644 --- a/CalPatRec/src/TZClusterFinder_module.cc +++ b/CalPatRec/src/TZClusterFinder_module.cc @@ -51,6 +51,7 @@ namespace mu2e { fhicl::Atom protonTCs {Name("protonTCs" ), Comment("Produce proton time clusters"), false}; fhicl::Atom useCCs {Name("useCCs" ), Comment("add CCs to TCs" ) }; fhicl::Atom recoverCCs {Name("recoverCCs" ), Comment("recover TCs using CCs" ) }; + fhicl::Atom requireCCs {Name("requireCCs" ), Comment("Require TCs with CCs" ), false }; fhicl::Atom chCollLabel {Name("chCollLabel" ), Comment("combo hit collection label" ) }; fhicl::Atom tcCollLabel {Name("tcCollLabel" ), Comment("time cluster coll label" ) }; fhicl::Atom ccCollLabel {Name("ccCollLabel" ), Comment("Calo Cluster coll label" ) }; @@ -87,6 +88,7 @@ namespace mu2e { bool _protonTCs; int _useCaloClusters; int _recoverCaloClusters; + bool _requireCaloClusters; std::optional> _timeoutService; std::optional _timeoutGuard; @@ -177,6 +179,7 @@ namespace mu2e { _protonTCs (config().protonTCs() ), _useCaloClusters (config().useCCs() ), _recoverCaloClusters (config().recoverCCs() ), + _requireCaloClusters (config().requireCCs() ), _timeoutService (std::nullopt ), _timeoutGuard (std::nullopt ), _chLabel (config().chCollLabel() ), @@ -312,6 +315,9 @@ namespace mu2e { //----------------------------------------------------------------------------- // put time cluster collection and protons counted into the event record //----------------------------------------------------------------------------- + if(_requireCaloClusters) { // remove time clusters without calo clusters from the output, keeping the counting/protons unchanged + std::erase_if(*tcColl, [](const TimeCluster& tc) { return !tc.hasCaloCluster();}); + } event.put(std::move(tcColl)); event.put(std::move(iiTC)); if(_protonTCs) event.put(std::move(protonTCColl), "proton"); diff --git a/CommonMC/src/SelectRecoMC_module.cc b/CommonMC/src/SelectRecoMC_module.cc index 01fa2504cf..37a258a80d 100644 --- a/CommonMC/src/SelectRecoMC_module.cc +++ b/CommonMC/src/SelectRecoMC_module.cc @@ -34,6 +34,7 @@ #include "Offline/RecoDataProducts/inc/ComboHit.hh" #include "Offline/RecoDataProducts/inc/CrvCoincidenceCluster.hh" #include "Offline/RecoDataProducts/inc/RecoCount.hh" +#include "Offline/RecoDataProducts/inc/TimeCluster.hh" // Utilities #include "Offline/ProditionsService/inc/ProditionsHandle.hh" #include "Offline/TrackerConditions/inc/StrawResponse.hh" @@ -60,25 +61,26 @@ namespace mu2e { using Name=fhicl::Name; using Comment=fhicl::Comment; struct Config { - fhicl::Atom debug { Name("debugLevel"), Comment("Debug Level"), 0}; - fhicl::Atom trkOnly { Name("TrkOnly"), Comment("Ignore cal and crv"), false}; - fhicl::Atom saveEnergySteps { Name("SaveEnergySteps"), Comment("Save all StepPoints that contributed energy to a StrawDigi"), false}; - fhicl::Atom saveUnused { Name("SaveUnusedDigiMCs"), Comment("Save StrawDigiMCs from particles used in any fit"), true}; - fhicl::Atom saveAllUnused { Name("SaveAllUnusedDigiMCs"), Comment("Save all StrawDigiMCs from particles used in any fit"), false}; - fhicl::Atom PP { Name("PrimaryParticle"), Comment("PrimaryParticle")}; - fhicl::Atom CCC { Name("CaloClusterCollection"), Comment("CaloClusterCollection")}; - fhicl::Sequence CrvCCCs { Name("CrvCoincidenceClusterCollections"),Comment("CrvCoincidenceClusterCollections")}; - fhicl::Atom SDC { Name("StrawDigiCollection"), Comment("StrawDigiCollection")}; - fhicl::Atom CHC { Name("ComboHitCollection"), Comment("ComboHitCollection for the original StrawHits (not Panel hits)")}; - fhicl::Atom CDC { Name("CaloDigiCollection"), Comment("CaloDigiCollection")}; - fhicl::Atom SDMCC { Name("StrawDigiMCCollection"), Comment("StrawDigiMCCollection")}; - fhicl::Atom CRVDC { Name("CrvDigiCollection"), Comment("CrvDigiCollection")}; - fhicl::Atom CRVDMCC { Name("CrvDigiMCCollection"), Comment("CrvDigiMCCollection")}; - fhicl::Atom PBTMC { Name("PBTMC"), Comment("ProtonBunchTimeMC")}; - fhicl::Atom EWM { Name("EventWindowMarker"), Comment("EventWindowMarker")}; - fhicl::Sequence KalSeeds { Name("KalSeedCollections"), Comment("KalSeedCollections")}; - fhicl::Atom VDSPC { Name("VDSPCollection"), Comment("Virtual Detector StepPointMC collection")}; - fhicl::Atom CCME { Name("CaloClusterMinE"), Comment("Minimum energy CaloCluster to save digis (MeV)")}; + fhicl::Atom debug { Name("debugLevel"), Comment("Debug Level"), 0}; + fhicl::Atom trkOnly { Name("TrkOnly"), Comment("Ignore cal and crv"), false}; + fhicl::Atom saveEnergySteps { Name("SaveEnergySteps"), Comment("Save all StepPoints that contributed energy to a StrawDigi"), false}; + fhicl::Atom saveUnused { Name("SaveUnusedDigiMCs"), Comment("Save StrawDigiMCs from particles used in any fit"), true}; + fhicl::Atom saveAllUnused { Name("SaveAllUnusedDigiMCs"), Comment("Save all StrawDigiMCs from particles used in any fit"), false}; + fhicl::Atom PP { Name("PrimaryParticle"), Comment("PrimaryParticle")}; + fhicl::Atom CCC { Name("CaloClusterCollection"), Comment("CaloClusterCollection")}; + fhicl::Sequence CrvCCCs { Name("CrvCoincidenceClusterCollections"),Comment("CrvCoincidenceClusterCollections")}; + fhicl::Atom SDC { Name("StrawDigiCollection"), Comment("StrawDigiCollection")}; + fhicl::Atom CHC { Name("ComboHitCollection"), Comment("ComboHitCollection for the original StrawHits (not Panel hits)")}; + fhicl::Atom CDC { Name("CaloDigiCollection"), Comment("CaloDigiCollection")}; + fhicl::Atom SDMCC { Name("StrawDigiMCCollection"), Comment("StrawDigiMCCollection")}; + fhicl::Atom CRVDC { Name("CrvDigiCollection"), Comment("CrvDigiCollection")}; + fhicl::Atom CRVDMCC { Name("CrvDigiMCCollection"), Comment("CrvDigiMCCollection")}; + fhicl::Atom PBTMC { Name("PBTMC"), Comment("ProtonBunchTimeMC")}; + fhicl::Atom EWM { Name("EventWindowMarker"), Comment("EventWindowMarker")}; + fhicl::Sequence KalSeeds { Name("KalSeedCollections"), Comment("KalSeedCollections")}; + fhicl::Sequence TimeClusters { Name("TimeClusterCollections"), Comment("TimeClusterCollections")}; + fhicl::Atom VDSPC { Name("VDSPCollection"), Comment("Virtual Detector StepPointMC collection")}; + fhicl::Atom CCME { Name("CaloClusterMinE"), Comment("Minimum energy CaloCluster to save digis (MeV)")}; }; using Parameters = art::EDProducer::Table; explicit SelectRecoMC(const Parameters& conf); @@ -103,7 +105,7 @@ namespace mu2e { bool _saveallenergy, _saveunused, _saveallunused; art::InputTag _pp, _ccc, _sdc, _chc, _cdc, _sdmcc, _crvdc, _crvdmcc, _pbtmc, _ewm, _vdspc; std::vector _kscs, _hscs; - std::vector _crvcccs; + std::vector _tccs, _crvcccs; double _ccme; // cache double _mbtime; // period of 1 microbunch @@ -132,6 +134,7 @@ namespace mu2e { _ewm(config().EWM()), _vdspc(config().VDSPC()), _kscs(config().KalSeeds()), + _tccs(config().TimeClusters()), _crvcccs(config().CrvCCCs()), _ccme(config().CCME()) { @@ -150,6 +153,9 @@ namespace mu2e { consumes(_crvdmcc); consumes(_pbtmc); consumes(_ewm); + for (const auto& i_tc : _tccs) { + consumes(i_tc); + } produces ("StrawDigiMap"); if (!_trkonly){ produces ("CrvDigiMap"); @@ -166,6 +172,7 @@ namespace mu2e { produces (); produces (); + if (_debug > 0) { std::cout << "Using KalSeed collections from "; @@ -435,6 +442,25 @@ namespace mu2e { std::set > simps; // add the MC primary SimParticles for(auto const& spp : pp.primarySimParticles()) simps.insert(spp); + + // save configured TimeCluster collections and preserve referenced straw/calo content + for (auto const& tcTag : _tccs) { + auto tcch = event.getHandle(tcTag); + if (!tcch.isValid()) continue; + auto const& tcc = *tcch; + + for (auto const& tc : tcc) { + for (auto const& hitIndex : tc.hits()) { + auto const ich = static_cast(hitIndex); + (void)chc.at(ich); // bounds check + sdindices.insert(static_cast(ich)); + } + if (tc.hasCaloCluster()) { + ccptrs.insert(tc.caloCluster()); + } + } + } + // loop over input KalFinalFit products for (auto const& ksc : _kscs) { // get all products from this diff --git a/CommonReco/src/SelectReco_module.cc b/CommonReco/src/SelectReco_module.cc index e484d81791..7fe77a307b 100644 --- a/CommonReco/src/SelectReco_module.cc +++ b/CommonReco/src/SelectReco_module.cc @@ -30,6 +30,7 @@ #include "Offline/RecoDataProducts/inc/KKLoopHelix.hh" #include "Offline/RecoDataProducts/inc/KKCentralHelix.hh" #include "Offline/RecoDataProducts/inc/KKLine.hh" +#include "Offline/RecoDataProducts/inc/TimeCluster.hh" // Utilities #include "Offline/ProditionsService/inc/ProditionsHandle.hh" #include "Offline/TrackerConditions/inc/StrawResponse.hh" @@ -53,17 +54,18 @@ namespace mu2e { using Name=fhicl::Name; using Comment=fhicl::Comment; struct Config { - fhicl::Atom debug { Name("debugLevel"), Comment("Debug Level"), 0}; - fhicl::Atom CCC { Name("CaloClusterCollection"), Comment("CaloClusterCollection")}; - fhicl::Atom CrvCCC { Name("CrvCoincidenceClusterCollection"),Comment("CrvCoincidenceClusterCollections")}; - fhicl::Atom SDC { Name("StrawDigiCollection"), Comment("StrawDigiCollection")}; - fhicl::Atom CHC { Name("ComboHitCollection"), Comment("ComboHitCollection for the original StrawHits (not Panel hits)")}; - fhicl::Atom CDC { Name("CaloDigiCollection"), Comment("CaloDigiCollection")}; - fhicl::Atom CRVDC { Name("CrvDigiCollection"), Comment("CrvDigiCollection")}; - fhicl::Sequence KalSeeds { Name("KalSeedCollections"), Comment("KalSeedCollections")}; + fhicl::Atom debug { Name("debugLevel"), Comment("Debug Level"), 0}; + fhicl::Atom CCC { Name("CaloClusterCollection"), Comment("CaloClusterCollection")}; + fhicl::Atom CrvCCC { Name("CrvCoincidenceClusterCollection"),Comment("CrvCoincidenceClusterCollections")}; + fhicl::Atom SDC { Name("StrawDigiCollection"), Comment("StrawDigiCollection")}; + fhicl::Atom CHC { Name("ComboHitCollection"), Comment("ComboHitCollection for the original StrawHits (not Panel hits)")}; + fhicl::Atom CDC { Name("CaloDigiCollection"), Comment("CaloDigiCollection")}; + fhicl::Atom CRVDC { Name("CrvDigiCollection"), Comment("CrvDigiCollection")}; + fhicl::Sequence KalSeeds { Name("KalSeedCollections"), Comment("KalSeedCollections")}; + fhicl::Sequence TimeClusters { Name("TimeClusterCollections"), Comment("TimeClusterCollections")}; // selection parameters - fhicl::Atom CCME { Name("CaloClusterMinE"), Comment("Minimum energy CaloCluster to save (MeV)")}; - fhicl::Atom saveNearby { Name("SaveNearbyDigis"), Comment("Save StrawDigis near the fit")}; + fhicl::Atom CCME { Name("CaloClusterMinE"), Comment("Minimum energy CaloCluster to save (MeV)")}; + fhicl::Atom saveNearby { Name("SaveNearbyDigis"), Comment("Save StrawDigis near the fit")}; }; using Parameters = art::EDProducer::Table; @@ -81,6 +83,7 @@ namespace mu2e { int debug_; art::InputTag ccct_, sdct_, chct_, cdct_, crvcct_, crvdct_; std::vector kscts_; + std::vector tccts_; double ccme_; bool saveunused_; }; @@ -95,6 +98,7 @@ namespace mu2e { crvcct_(config().CrvCCC()), crvdct_(config().CRVDC()), kscts_(config().KalSeeds()), + tccts_(config().TimeClusters()), ccme_(config().CCME()), saveunused_(config().saveNearby()) { @@ -106,6 +110,10 @@ namespace mu2e { consumes(crvdct_); consumes(crvcct_); consumesMany(); + for (auto const& tcct : tccts_) { + consumes(tcct); + } + produces ("StrawDigiMap"); produces ("CrvDigiMap"); produces (); @@ -115,10 +123,15 @@ namespace mu2e { produces (); produces (); produces (); + for (auto const& tcct : tccts_) { + produces(tcct.label()); + } if (debug_ > 0) { std::cout << "Using KalSeed collections from "; for (auto const& kff : kscts_) std::cout << kff << " " << std::endl; + std::cout << "Using TimeCluster collections from "; + for (auto const& tcc : tccts_) std::cout << tcc << " " << std::endl; } } @@ -134,7 +147,7 @@ namespace mu2e { } void SelectReco::fillTrk( art::Event& event, std::set >& ccptrs, RecoCount& nrec) { - // Tracker-reated data products + // Tracker-reated data products auto sdch = event.getValidHandle(sdct_); auto const& sdc = *sdch; auto sdadcch = event.getValidHandle(sdct_); @@ -144,10 +157,35 @@ namespace mu2e { // create products related to the reconstruction output std::unique_ptr ssdc(new StrawDigiCollection); std::unique_ptr ssdadcc(new StrawDigiADCWaveformCollection); + std::vector> stccs; + stccs.reserve(tccts_.size()); // index maps between original collections and pruned collections std::unique_ptr sdim(new IndexMap); // straw digi indices that are referenced by the tracks, or are 'close to' the track SDIS sdindices; + + // save configured TimeCluster collections and keep referenced content + for (auto const& tcct : tccts_) { + std::unique_ptr stcc(new TimeClusterCollection); + auto tcch = event.getHandle(tcct); + if (tcch.isValid()) { + auto const& tcc = *tcch; + stcc->reserve(tcc.size()); + for (auto const& tc : tcc) { + for (auto const& hitIndex : tc.hits()) { + auto const ich = static_cast(hitIndex); + (void)chc.at(ich); // bounds check vs input ComboHitCollection + sdindices.insert(static_cast(ich)); + } + if (tc.hasCaloCluster()) { + ccptrs.insert(tc.caloCluster()); + } + stcc->push_back(tc); + } + } + stccs.push_back(std::move(stcc)); + } + // loop over input KalSeeds for (auto const& ksct : kscts_){ auto ksch = event.getHandle(ksct); @@ -179,6 +217,9 @@ namespace mu2e { event.put(std::move(sdim),"StrawDigiMap"); event.put(std::move(ssdc)); event.put(std::move(ssdadcc)); + for (size_t i = 0; i < tccts_.size(); ++i) { + event.put(std::move(stccs.at(i)), tccts_.at(i).label()); + } } void SelectReco::fillStrawHitCounts(ComboHitCollection const& chc, RecoCount& nrec) { diff --git a/CosmicReco/src/LineFinder_module.cc b/CosmicReco/src/LineFinder_module.cc index 9a7a865246..b79e355493 100644 --- a/CosmicReco/src/LineFinder_module.cc +++ b/CosmicReco/src/LineFinder_module.cc @@ -16,12 +16,14 @@ #include "Offline/GeometryService/inc/GeomHandle.hh" #include "Offline/GeometryService/inc/DetectorSystem.hh" +#include "Offline/CalorimeterGeom/inc/Calorimeter.hh" #include "Offline/TrackerGeom/inc/Tracker.hh" #include "Offline/Mu2eUtilities/inc/TwoLinePCA.hh" #include "Offline/RecoDataProducts/inc/ComboHit.hh" #include "Offline/RecoDataProducts/inc/TimeCluster.hh" #include "Offline/RecoDataProducts/inc/CosmicTrackSeed.hh" +#include "Offline/RecoDataProducts/inc/CaloCluster.hh" #include "CLHEP/Units/PhysicalConstants.h" #include "CLHEP/Matrix/Vector.h" @@ -56,11 +58,14 @@ namespace mu2e{ fhicl::Atom nmax{Name("MaxPairs"), Comment("Max pairs to try")}; fhicl::Atom chToken{Name("ComboHitCollection"),Comment("tag for straw hit collection")}; fhicl::Atom tcToken{Name("TimeClusterCollection"),Comment("tag for time cluster collection")}; + fhicl::Atom ccToken{Name("CaloClusterCollection"),Comment("tag for calo cluster collection"), ""}; + fhicl::Atom addCaloClusters{Name("AddCaloClusters"), Comment("Whether to connect calo clusters to the line seeds"), false}; }; typedef art::EDProducer::Table Parameters; explicit LineFinder(const Parameters& conf); virtual ~LineFinder(){}; virtual void produce(art::Event& event ) override; + virtual void beginRun(art::Run& run) override; private: @@ -76,11 +81,17 @@ namespace mu2e{ unsigned _nmax; art::InputTag _chToken; art::InputTag _tcToken; + art::InputTag _ccToken; + bool _addCaloClusters; ProditionsHandle _alignedTracker_h; + const CaloClusterCollection* _cccol; + art::Handle _cccolH; + const Calorimeter* _calorimeter; int findLine(const ComboHitCollection& shC, std::vector const& shiv, art::Event const& event, CosmicTrackSeed& tseed); - }; + void addClusterToSeed(CosmicTrackSeed& tseed); +}; LineFinder::LineFinder(const Parameters& conf) : @@ -93,15 +104,26 @@ namespace mu2e{ _Ntsteps (conf().ntsteps()), _stepSize (conf().stepsize()), _nmax (conf().nmax()), - _chToken (conf().chToken()), - _tcToken (conf().tcToken()) + _chToken (conf().chToken()), + _tcToken (conf().tcToken()), + _ccToken (conf().ccToken()), + _addCaloClusters (conf().addCaloClusters()), + _cccol(nullptr), + _calorimeter(nullptr) { consumes(_chToken); consumes(_tcToken); + if(_addCaloClusters) { + consumes(_ccToken); + } produces(); } +void LineFinder::beginRun(art::Run& run) { + _calorimeter = GeomHandle().get(); +} + void LineFinder::produce(art::Event& event ) { auto const& chH = event.getValidHandle(_chToken); @@ -109,9 +131,18 @@ void LineFinder::produce(art::Event& event ) { auto const& tcH = event.getValidHandle(_tcToken); const TimeClusterCollection& tccol(*tcH); + if(_addCaloClusters) { + event.getByLabel(_ccToken, _cccolH); + if(!_cccolH.isValid()) { + throw cet::exception("RECO") << "[LineFinder::produce] Could not retrieve calo cluster collection with tag " << _ccToken.encode(); + } + _cccol = _cccolH.product(); + } std::unique_ptr seed_col(new CosmicTrackSeedCollection()); + if(_diag > 2) std::cout << "[LineFinder:: " << __func__ << "] Event " << event.id() + << " has " << tccol.size() << " time clusters\n"; for (size_t index=0;index< tccol.size();++index) { const auto& tclust = tccol[index]; @@ -126,13 +157,22 @@ void LineFinder::produce(art::Event& event ) { tseed._timeCluster = art::Ptr(tcH,index); tseed._track.converged = true; + if(_diag > 2) std::cout << " Time cluster " << index << " has " << tclust.hits().size() << " hits" + << " and " << chc.size() << " straw hits\n"; int seedSize = findLine(chc, shiv, event, tseed); if (_diag > 1) - std::cout << "LineFinder: seedSize = " << seedSize << std::endl; + std::cout << "[LineFinder::" << __func__ <<"] Line seed N(hits) = " << seedSize << std::endl; if (seedSize >= _minPeak){ if (_diag > 0) - std::cout << "LineFinder: found line (" << seedSize << ")" << " " << tseed._track.FitParams.T0 << " " << tseed._track.FitParams.A0 << " " << tseed._track.FitParams.B0 << " " << tseed._track.FitParams.A1 << " " << tseed._track.FitParams.B1 << std::endl; + std::cout << "LineFinder: found line (N(hits) = " << seedSize << ")" + << " " << tseed._track.FitParams.T0 + << " " << tseed._track.FitParams.A0 + << " " << tseed._track.FitParams.B0 + << " " << tseed._track.FitParams.A1 + << " " << tseed._track.FitParams.B1 + << std::endl; + if(_addCaloClusters) addClusterToSeed(tseed); tseed._status.merge(TrkFitFlag::Straight); tseed._status.merge(TrkFitFlag::hitsOK); tseed._status.merge(TrkFitFlag::helixOK); @@ -271,6 +311,51 @@ int LineFinder::findLine(const ComboHitCollection& shC, std::vectorsize(); ++i_cl) { + const CaloCluster& cl = _cccol->at(i_cl); + const CLHEP::Hep3Vector cl_pos = _calorimeter->geomUtil().mu2eToTracker(_calorimeter->geomUtil().diskToMu2e(cl.diskID(), cl.cog3Vector())); + const double cl_time = cl.time(); + // Get the expected time at the calo disk based on the line seed + const double speed = CLHEP::c_light; // assume speed of the particle is close to the speed of light + const double dr = std::abs((cl_pos - seedIntCLHEP).dot(seedDirCLHEP.unit())); + const double t_flight = dr / speed; + // check either trajectory + const double t0_diff_1 = std::abs((t0 + t_flight) - cl_time); + const double t0_diff_2 = std::abs((t0 - t_flight) - cl_time); + const double t0_diff = std::min(t0_diff_1, t0_diff_2); + if(t0_diff < min_t0_diff && t0_diff < t0_match_window) { + min_t0_diff = t0_diff; + matchedCluster = &cl; + matchedClusterIndex = i_cl; + } + } + + if(matchedCluster) { + tseed._caloCluster = art::Ptr(_cccolH, matchedClusterIndex); + if(_diag > 0) { + std::cout << "LineFinder: matched calo cluster with energy " << matchedCluster->energyDep() << " MeV" + << " time " << matchedCluster->time() << " ns to line seed with t0 = " + << tseed._t0._t0 << " ns" << std::endl; + } + } +} + }//end mu2e namespace using mu2e::LineFinder; DEFINE_ART_MODULE(LineFinder) diff --git a/RecoDataProducts/inc/CosmicTrackSeed.hh b/RecoDataProducts/inc/CosmicTrackSeed.hh index bde53c0f60..d9b28ecb3d 100644 --- a/RecoDataProducts/inc/CosmicTrackSeed.hh +++ b/RecoDataProducts/inc/CosmicTrackSeed.hh @@ -6,6 +6,7 @@ #include "Offline/RecoDataProducts/inc/ComboHit.hh" #include "Offline/RecoDataProducts/inc/TrkStrawHitSeed.hh" #include "Offline/RecoDataProducts/inc/CosmicTrack.hh" +#include "Offline/RecoDataProducts/inc/CaloCluster.hh" #include "Offline/RecoDataProducts/inc/TimeCluster.hh" #include "Offline/RecoDataProducts/inc/TrkFitFlag.hh" #include "canvas/Persistency/Common/Ptr.h" @@ -21,13 +22,16 @@ namespace mu2e { CosmicTrack const& track() const { return _track; } TrkFitFlag const& status() const { return _status; } art::Ptr const& timeCluster() const { return _timeCluster; } + art::Ptr const& caloCluster() const { return _caloCluster; } ComboHitCollection const& hits() const { return _straw_chits; } + bool hasCaloCluster() const { return _caloCluster.isNonnull(); } TrkT0 _t0; // t0 for this track CosmicTrack _track; // Cosmic track created from these hits TrkFitFlag _status; // status of processes used to create this seed // helixOK: ???, helixConverged: ???, circleInit: ???, Straight: ???, hitsOK: ??? art::Ptr _timeCluster; // associated time cluster + art::Ptr _caloCluster; // associated calo cluster ComboHitCollection _straw_chits; // get the straw level hits and store here (need to find panel hits first) // For future use only