Add optional proton time cluster output, update printouts

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

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);
     }
   }

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);

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) {

CalPatRec/src/TZClusterFinder_module.cc

@@ -48,6 +48,7 @@ namespace mu2e {
       fhicl::Atom<int>               debugLevel       {Name("debugLevel"       ), Comment("turn on/off debug"           ) };
       fhicl::Atom<int>               printFrequency   {Name("printFrequency"   ), Comment("print frequency"             ) };
       fhicl::Atom<int>               runDisplay       {Name("runDisplay"       ), Comment("will plot t vs z"            ) };
+      fhicl::Atom<bool>              protonTCs        {Name("protonTCs"        ), Comment("Produce proton time clusters"), false};
       fhicl::Atom<int>               useCCs           {Name("useCCs"           ), Comment("add CCs to TCs"              ) };
       fhicl::Atom<int>               recoverCCs       {Name("recoverCCs"       ), Comment("recover TCs using CCs"       ) };
       fhicl::Atom<art::InputTag>     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<art::ServiceHandle<TimeoutWatchdogService>> _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<CaloClusterCollection>(_ccLabel);
       produces<TimeClusterCollection>();
       produces<IntensityInfoTimeCluster>();
+      if(_protonTCs) produces<TimeClusterCollection>("proton");
 
 
       if (_runDisplay == 1) { _c1 = new TCanvas("_c1", "t vs. z", 900, 900); }
@@ -290,9 +294,11 @@ namespace mu2e {
 
     std::unique_ptr<IntensityInfoTimeCluster> iiTC(new IntensityInfoTimeCluster);
     std::unique_ptr<TimeClusterCollection>    tcColl(new TimeClusterCollection);
+    std::unique_ptr<TimeClusterCollection>    protonTCColl = (_protonTCs) ? std::make_unique<TimeClusterCollection>() : 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<mu2e::CaloCluster>(_ccHandle, caloIdx);
+        tc._nsh = _f.chunks[i].nStrawHits;
+        tcs->push_back(std::move(tc));
+      }
     }
 
     outIITC.setNProtonTCs(nProtons);