ActivityTrigger_module.cc

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////
// Class:       ActivityTrigger
// Module Type: filter
// File:        ActivityTrigger_module.cc
//
// Generated at Wed Sep 12 16:15:01 2012 by Nova DDT using artmod
// from art v1_01_01.
////////////////////////////////////////////////////////////////////////

#include "art/Framework/Core/EDFilter.h"
#include "art/Framework/Core/ModuleMacros.h"
#include "canvas/Persistency/Common/Assns.h"
#include "art/Framework/Principal/Handle.h"
#include "art/Framework/Principal/Event.h"

#include "DDTBaseDataProducts/BaseProducts.h"
#include "DDTBaseDataProducts/HitList.h"
#include "DDTBaseDataProducts/DAQHit.h"
#include "DDTBaseDataProducts/TriggerDecision.h"
#include "DDTBaseDataProducts/GroupedHitList.h"
#include "DDTUtilities/AssociationUtil.h"

#include "DAQDataFormats/DAQDataFormats.h"

#include <string>
#include <algorithm>
#include <unordered_set>

namespace novaddt {
  class ActivityTrigger;
}

class novaddt::ActivityTrigger : public art::EDFilter {

public:
  explicit ActivityTrigger(fhicl::ParameterSet const & p);
  virtual ~ActivityTrigger();

  virtual bool filter(art::Event & e);

private:

  std::string  _hitsModuleLabel;   ///< label of module making the HitList
  std::string  _hitsInstanceLabel; ///< instance label making the HitList
  unsigned int _minHits;           ///< Minimum number of hits in the slice
  unsigned int _minPlanesView;     ///< Minimum planes (unique) on each view
  unsigned int _minPlanes;         ///< Minimum planes (unique) the hits are spread in any view
  unsigned int _maxNearWindow;     ///< Maximum window where to look for near planes
  unsigned int _minNearPlanes;     ///< Minimum number of near planes
  unsigned int _prescale;          ///< holds prescaling variable from fcl file
  
  unsigned int _triggerCounts = 0;
  unsigned int _evtCounts = 0;
};

//---------------------------------------------------------------------    
novaddt::ActivityTrigger::ActivityTrigger(fhicl::ParameterSet const & p)
  : _hitsModuleLabel  (p.get<std::string >("HitModuleLabel"  ))
  , _hitsInstanceLabel(p.get<std::string >("HitInstanceLabel"))
  , _minHits          (p.get<unsigned int>("MinHits"         ))
  , _minPlanesView    (p.get<unsigned int>("MinPlanesView"   ))
  , _minPlanes        (p.get<unsigned int>("MinPlanes"       ))
  , _maxNearWindow    (p.get<unsigned int>("MaxNearWindow"   ))
  , _minNearPlanes    (p.get<unsigned int>("MinNearPlanes"   ))
  , _prescale         (p.get<unsigned int>("prescale"        ))
{
  produces< std::vector<novaddt::TriggerDecision> >();
  //produces< art::Assns<novaddt::TriggerDecision, novaddt::HitList> >();
}

//---------------------------------------------------------------------    
novaddt::ActivityTrigger::~ActivityTrigger()
{
  // Clean up dynamic memory and other resources here.
}

//---------------------------------------------------------------------    
bool novaddt::ActivityTrigger::filter(art::Event & e)
{
  //std::cout << "--- ActivityTrigger filter\n";
  ++_evtCounts;

  // Retrieve the (grouped) hit list for this event
  art::Handle<novaddt::GroupedHitList<novaddt::TDC> > hits;
  e.getByLabel(_hitsModuleLabel, _hitsInstanceLabel, hits);

  // Allocate the vector for the trigger decisions that we may make (and write to the event)
  std::unique_ptr<std::vector<novaddt::TriggerDecision>> td(new std::vector<novaddt::TriggerDecision>);

  /**************************************************************************
   * Activity Trigger Algorithm                                             *
   *                                                                        *
   * Essentially all we do is count the number of hits we have in slice     *
   * if it is over a given threshold then we write it out.                  *
   *                                                                        *
   * Overlaps with the beam (numi) trigger are handled via an inhibit at    *
   * the global trigger level                                               *
   *                                                                        *
   * However there are a couple of additional cuts we put on the data,      *
   * first in addition to having at least N hits we require that they       *
   * extend over at least Mx planes in the X view, and analogously with     *
   * the Y view. We also require the total number of planes to be larger    *
   * than P. Finally, we perform a continuity check, asking for K planes    *
   * with hits whithin a W window.                                          *
   *                                                                        *
   * This allows the trigger to be tuned in the following manners           *
   * Increase: N = more hits (more activity, larger showers)                *
   * Increase: Mx and P = more "horizontal" (z-aligned) activity            *
   * Decrease M or P = more "vertical" (x/y-aligned activity)               *
   * Increase K: longer continuous tracks                                   *
   * Decrease W: less gaps are allowed. If W<K or any of them is set to 0   *
   *  no continuity check is performed                                      *
   **************************************************************************/

  std::unordered_set<int> XplaneList;
  std::unordered_set<int> YplaneList;
  std::unordered_set<int> planeList;

  // Loop over the slices in the event
  for(unsigned int i=0; i<hits->n_groups(); ++i){

    // Get the hits for the current slice
    novaddt::HitList theSliceHits = hits->get(i);

    XplaneList.clear();
    YplaneList.clear();
    planeList.clear();
    bool fContinuity = false;
    // Safety checks in case we mess up the fickle file
    if(_maxNearWindow < _minNearPlanes) fContinuity = true; 
    if(!(_maxNearWindow * _minNearPlanes)) fContinuity = true;

    // Test to see how many hits are in the slice
    if (theSliceHits.size() >= _minHits){
      // Fill the unique planes
      for(unsigned int j=0; j < theSliceHits.size(); ++j){
        auto thePlane = theSliceHits.at(j).Plane().val;
        auto theView = theSliceHits.at(j).View().val;
        // Push the plane into list of unique planes in the proper view
        if(theView == daqchannelmap::X_VIEW) XplaneList.insert(thePlane);
        else if(theView == daqchannelmap::Y_VIEW) YplaneList.insert(thePlane);
        planeList.insert(thePlane);
      } // endfor j
      
      if((XplaneList.size() + YplaneList.size()) >= _minPlanes){
        // If there are enough planes, check by view
        if(XplaneList.size() >= _minPlanesView){
          if(YplaneList.size() >= _minPlanesView){
            /**********************************************************
             * Check for continuity. This means a minimum number of    *
             * planes _minNearPlanes within a window of _maxNearPlanes *
             * (e.g., 4 out of 5 planes, a la MINOS)                   *
             **********************************************************/
            for(auto iPlane=planeList.begin(); iPlane!=planeList.end(); ++iPlane)
              {
                if(!fContinuity)
                  {
                    unsigned int fNearPlanes = 0;
                    for(unsigned int j=1; j<=_maxNearWindow; ++j)
                      fNearPlanes += planeList.count(*iPlane+j);
                    
                    if(fNearPlanes >= _minNearPlanes)
                      fContinuity = true;
                  }
              } // end of continuity check
            
            if(fContinuity)
              {
                /******************************************************
                 * EVENTS THAT MAKE IT HERE SHOULD GENERATE A TRIGGER *
                 ******************************************************/
                // Internal prescale 
                ++_triggerCounts;
                if (_triggerCounts % _prescale == _prescale - 1){
                  // Store a trigger decision for this slice
                  td->emplace_back(theSliceHits.front().TDC().val, 
                                   theSliceHits.back().TDC().val - theSliceHits.front().TDC().val, 
                                   daqdataformats::trigID::TRIG_ID_DATA_ACTIVITY_1,
                                   _prescale);
                } // end of issue trigger decision
              } // end of cointinuity check       
          } // endif gt minYPlanes
        } // endif gt minXPlanes
      } // endif gt minPlanes
    } // end of check on hit list size
  } // end loop on hit groups

  bool goodTrigger = (td->size() > 0);

  e.put(std::move(td));
  //e.put(std::move(assn));

  //std::cout << "Frequency: " << _triggerCounts*1./(_evtCounts*5.e-3) << " pm " // Fraction of passing per event. One event lasts 5 ms 
  //          << sqrt(_triggerCounts*(_evtCounts-_triggerCounts)*1./_evtCounts)/(_evtCounts*5.e-3) << " Hz" << std::endl; // Binomial Error
  
  return goodTrigger;
}

DEFINE_ART_MODULE(novaddt::ActivityTrigger)