diff --git a/CalPatRec/inc/TZClusterFinder_types.hh b/CalPatRec/inc/TZClusterFinder_types.hh index b76d726be3..bbffac755d 100644 --- a/CalPatRec/inc/TZClusterFinder_types.hh +++ b/CalPatRec/inc/TZClusterFinder_types.hh @@ -71,6 +71,7 @@ namespace mu2e { const CaloClusterCollection* _ccColl; TimeClusterCollection* _tcColl; // 'tcColl': time cluster collection IntensityInfoTimeCluster* _iiTC; + TimeClusterCollection* _protonTCColl = nullptr; // proton time clusters, optionally filled int _nTZClusters; // diagnostic data members used in TZ tool diff --git a/CalPatRec/src/AgnosticHelixFinderDiagDisplay.cc b/CalPatRec/src/AgnosticHelixFinderDiagDisplay.cc index c6432a0d75..956cb745c7 100644 --- a/CalPatRec/src/AgnosticHelixFinderDiagDisplay.cc +++ b/CalPatRec/src/AgnosticHelixFinderDiagDisplay.cc @@ -774,9 +774,13 @@ void AgnosticHelixFinderDiag::plotHelixStageXY(int stage) { } else { // All helices const bool use_tc = stage == kTimeCluster; if(use_tc) title = "Helices (time cluster)"; + if(_debugLevel > 1) printf(" %12s: Plotting %zu helices, stage = %i\n", + __func__, _data->helixSeedData.size(), stage); plotHitsXY(use_tc); // draw (relevant) hits plotMC(kHelix); for(size_t index = 0; index < _data->helixSeedData.size(); ++index) { + if(_debugLevel > 1) printf(" %12s: Plotting the helix XY circle %zu, stage = %i\n", + __func__, index, stage); plotHelixXY(_data->helixSeedData.at(index).seed, index); } } diff --git a/CalPatRec/src/AgnosticHelixFinderDiag_tool.cc b/CalPatRec/src/AgnosticHelixFinderDiag_tool.cc index 6afbca9aae..d9c2cadb0e 100644 --- a/CalPatRec/src/AgnosticHelixFinderDiag_tool.cc +++ b/CalPatRec/src/AgnosticHelixFinderDiag_tool.cc @@ -515,7 +515,7 @@ int AgnosticHelixFinderDiag::fillHistograms(void* Data, int Mode) { if(_data->hseed) { // only display at this level, where per-helix wasn't displayed but accepted helices are each plotted if(_display && _data->diagLevel == 3) { - std::cout << "Helix stage: Loop condition " << ConditionName(_data->loopCondition) + std::cout << "Final stage: Loop condition " << ConditionName(_data->loopCondition) << std::endl; plotTripletStage(); plotHelixStageXY(kHelix); @@ -535,6 +535,8 @@ int AgnosticHelixFinderDiag::fillHistograms(void* Data, int Mode) { if(Mode == DIAG::kTimeCluster) { if(_display && _data->diagLevel > 1) { + std::cout << "TimeCluster stage:" + << std::endl; // plotTripletStage(true); plotCircleStage(true); plotHelixStageXY(kTimeCluster); @@ -550,6 +552,8 @@ int AgnosticHelixFinderDiag::fillHistograms(void* Data, int Mode) { if(Mode == DIAG::kEnd) { if(_display) { + std::cout << "End stage:" + << std::endl; // plotTripletStage(true); plotCircleStage(true); plotHelixStageXY(kEnd); diff --git a/CalPatRec/src/AgnosticHelixFinder_module.cc b/CalPatRec/src/AgnosticHelixFinder_module.cc index d4c92b582e..195329d7a4 100644 --- a/CalPatRec/src/AgnosticHelixFinder_module.cc +++ b/CalPatRec/src/AgnosticHelixFinder_module.cc @@ -500,6 +500,8 @@ namespace mu2e { const int nHelicesInitial = _diagInfo.nHelices; // Only valid if diagLevel > 0, for diagnostic tracking tcHitsFill(i); // Initialize the list of hits in the time cluster while(findHelix(i, *hsColl) && _findMultipleHelices); // Exit the search if no helix is found or after finding a helix if not configured for multi-helix reco + if(_diagLevel > 1) + printf(" Time cluster %zu has %i helices\n", i, _diagInfo.nHelices - nHelicesInitial); if (_diagLevel > 0) { tcInfo timeClusterInfo; @@ -515,6 +517,8 @@ namespace mu2e { printf("[AgnosticHelixFinder::%s] %i:%i:%i Event data not found!\n", __func__, event.run(), event.subRun(), event.event()); } if(_doTiming) _watch->StopTime("per-time cluster"); + if(_diagLevel > 1) + printf(" Found %i helices total (%zu in collection)\n", _diagInfo.nHelices, hsColl->size()); // fill necessary data members for diagnostic tool if (_diagLevel > 0) { diff --git a/CalPatRec/src/TZClusterFinder_module.cc b/CalPatRec/src/TZClusterFinder_module.cc index 39a7303c33..fe4b465cb6 100644 --- a/CalPatRec/src/TZClusterFinder_module.cc +++ b/CalPatRec/src/TZClusterFinder_module.cc @@ -48,6 +48,7 @@ namespace mu2e { fhicl::Atom debugLevel {Name("debugLevel" ), Comment("turn on/off debug" ) }; fhicl::Atom printFrequency {Name("printFrequency" ), Comment("print frequency" ) }; fhicl::Atom runDisplay {Name("runDisplay" ), Comment("will plot t vs z" ) }; + 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 chCollLabel {Name("chCollLabel" ), Comment("combo hit collection label" ) }; @@ -83,6 +84,7 @@ namespace mu2e { int _debugLevel; int _printfreq; int _runDisplay; + bool _protonTCs; int _useCaloClusters; int _recoverCaloClusters; std::optional> _timeoutService; @@ -172,10 +174,11 @@ namespace mu2e { _debugLevel (config().debugLevel() ), _printfreq (config().printFrequency() ), _runDisplay (config().runDisplay() ), + _protonTCs (config().protonTCs() ), _useCaloClusters (config().useCCs() ), _recoverCaloClusters (config().recoverCCs() ), - _timeoutService (std::nullopt ), - _timeoutGuard (std::nullopt ), + _timeoutService (std::nullopt ), + _timeoutGuard (std::nullopt ), _chLabel (config().chCollLabel() ), _tcLabel (config().tcCollLabel() ), _ccLabel (config().ccCollLabel() ), @@ -207,6 +210,7 @@ namespace mu2e { consumes(_ccLabel); produces(); produces(); + if(_protonTCs) produces("proton"); if (_runDisplay == 1) { _c1 = new TCanvas("_c1", "t vs. z", 900, 900); } @@ -290,9 +294,11 @@ namespace mu2e { std::unique_ptr iiTC(new IntensityInfoTimeCluster); std::unique_ptr tcColl(new TimeClusterCollection); + std::unique_ptr protonTCColl = (_protonTCs) ? std::make_unique() : nullptr; _data._tcColl = tcColl.get(); _data._iiTC = iiTC.get(); + _data._protonTCColl = (_protonTCs) ? protonTCColl.get() : nullptr; bool ok = findData(event); @@ -306,6 +312,7 @@ namespace mu2e { //----------------------------------------------------------------------------- event.put(std::move(tcColl)); event.put(std::move(iiTC)); + if(_protonTCs) event.put(std::move(protonTCColl), "proton"); _timeoutGuard.reset(); @@ -756,10 +763,24 @@ namespace mu2e { unsigned short nProtons = 0; + auto tcs = (_protonTCs) ? _data._protonTCColl : nullptr; for (size_t i=0; i<_f.chunks.size(); i++) { if (_f.chunks[i].nrgSelection == 1) {continue;} if (_f.chunks[i].nStrawHits < _clusterThresh) {continue;} nProtons++; + + // If proton time clusters are requested, add them to the output collection + if(_protonTCs) { + TimeCluster tc; + for (size_t j=0; j<_f.chunks[i].hIndices.size(); j++) { + tc._strawHitIdxs.push_back(StrawHitIndex(_f.chunks[i].hIndices[j])); + } + tc._t0 = TrkT0(_f.chunks[i].fitter.y0(), 0.); + const int caloIdx = _f.chunks[i].caloIndex; + if (_useCaloClusters == 1 && caloIdx != -1) tc._caloCluster = art::Ptr(_ccHandle, caloIdx); + tc._nsh = _f.chunks[i].nStrawHits; + tcs->push_back(std::move(tc)); + } } outIITC.setNProtonTCs(nProtons);