JTriggerMXShower.hh

Go to the documentation of this file.

#ifndef __JTRIGGER__JTRIGGERMXSHOWER__
#define __JTRIGGER__JTRIGGERMXSHOWER__
 
#include <vector>
#include <set>
 
#include "JTrigger/JAlgorithm.hh"
#include "JTrigger/JMatch3G.hh"
#include "JTrigger/JParametersHelper.hh"
#include "JTrigger/JTriggerInterface.hh"
#include "JTrigger/JTimeslice.hh"
#include "JTrigger/JTimesliceClone.hh"
#include "JTrigger/JTriggerOutput.hh"
#include "JTrigger/JTriggerMXShower_t.hh"
#include "JDetector/JDetector.hh"
#include "JDetector/JModuleMapper.hh"
#include "JMath/JMathToolkit.hh"
#include "JPhysics/JConstants.hh"
#include "JGeometry3D/JVector3D.hh"
 
 
/**
 * \author mdejong
 */
 
namespace JTRIGGER {}
namespace JPP { using namespace JTRIGGER; }
 
namespace JTRIGGER {
 
  using JDETECTOR::JModule;
  using JDETECTOR::JModuleMapper;
  using JDETECTOR::JDetector;
  using JDETECTOR::JMaximalDistance;
  using JGEOMETRY3D::JVector3D;
 
 
  /**
   * Module attributes.
   */
  struct JAttributes {
    double Tmin_ns;         //!< minimal time difference [ns]
    double Tmax_ns;         //!< maximal time difference [ns]
  };
 
 
  /**
   * Set module attributes.
   *
   * \param  first          first  module
   * \param  second         second module
   * \param  attributes     module attributes
   */
  inline void setAttributes(const JModule& first,
                            const JModule& second,
                            JAttributes&   attributes)
  {
    using namespace JPP;
 
    typedef JTriggerMXShower_t::JParameters  JParameters_t;
 
    // get trigger parameters
 
    const JParameters_t& parameters = JParametersHelper<JParameters_t>::getParameters();
 
    const double d = getDistance(first.getPosition(), second.getPosition());
 
    double t1 = 0.0;
 
    if      (d <= 0.5 * parameters.DMax_m)
      t1 =  d                      * getIndexOfRefraction() * getInverseSpeedOfLight()  +  parameters.TMaxExtra_ns;
    else if (d <= parameters.DMax_m)
      t1 = (parameters.DMax_m - d) * getIndexOfRefraction() * getInverseSpeedOfLight()  +  parameters.TMaxExtra_ns;
 
    attributes.Tmin_ns = -t1;
    attributes.Tmax_ns = +t1;
  }
 
 
  /**
   * Shower trigger.
   */
  class JTriggerMXShower : 
    public JTriggerMXShower_t,
    public JTriggerInterface,
    public std::vector<JVector3D>
  {
  public:
 
    /**
     * Type definition of module mapper.
     */
    typedef JModuleMapper<JAttributes>   JModuleMapper_t;
 
 
    /**
     * Constructor.
     *
     * \param  input          trigger parameters
     * \param  detector       detector;
     */
    JTriggerMXShower(const JParameters& input,
                     const JDetector&   detector) :
      parameters(input),
      mapper(detector, JMaximalDistance(parameters.DMax_m)),
      match(parameters.DMax_m, parameters.TMaxExtra_ns)
    {
      using namespace JPP;
 
      TMaxEvent_ns = parameters.DMax_m * getIndexOfRefraction() * getInverseSpeedOfLight()  +  parameters.TMaxExtra_ns;
 
      if (parameters.numberOfHits == 0) {
        parameters.numberOfHits = 1;
      }
 
      if (parameters.numberOfHits == 1) {
        parameters.factoryLimit = 1;
      }
 
      if (parameters.numberOfModules < 1) {
        parameters.numberOfModules = 1;
      }
 
      // vertices
 
      if (parameters.RMax_m > 0.0 &&
          parameters.Xv_m   > 0.0) {
 
        this->push_back(JVector3D(0.0,0.0,0.0));
 
        for (double x = 0.5*parameters.Xv_m; x < parameters.RMax_m + 0.5*parameters.Xv_m; x += parameters.Xv_m) {
          for (double y = 0.5*parameters.Xv_m; y < parameters.RMax_m + 0.5*parameters.Xv_m; y += parameters.Xv_m) {
            for (double z = 0.5*parameters.Xv_m; z < parameters.RMax_m + 0.5*parameters.Xv_m; z += parameters.Xv_m) {
              if (x*x + y*y + z*z <= parameters.RMax_m*parameters.RMax_m) {
                this->push_back(JVector3D(-x,-y,-z));
                this->push_back(JVector3D(-x,-y,+z));
                this->push_back(JVector3D(-x,+y,-z));
                this->push_back(JVector3D(-x,+y,+z));
                this->push_back(JVector3D(+x,-y,-z));
                this->push_back(JVector3D(+x,-y,+z));
                this->push_back(JVector3D(+x,+y,-z));
                this->push_back(JVector3D(+x,+y,+z));
              }
            }
          }
        }
      }
    }
 
 
    /**
     * Process trigger.
     *
     * \param  inputL1        input L1 data
     * \param  inputL0        input L0 data
     * \param  out            output data
     */
    template<class JElement_t>
    void operator()(const JTriggerInput&         inputL1,
                    const JTimeslice<JElement_t> inputL0,
                    std::back_insert_iterator<JTriggerOutput> out) const
    {
      typedef JTimesliceClone<JElement_t>                JTimeslice_t;
      typedef typename JTimeslice_t::value_type          JFrame_t;
      typedef typename JModuleMapper_t::container_type   container_type;
 
      if (parameters.enabled) {
 
        JTimeslice_t clone(inputL0, mapper);
 
        for (JTriggerInput::const_iterator root = inputL1.begin(); root != inputL1.end(); ++root) {
 
          buffer.clear();
 
          const container_type& zip = mapper.getList(root->getModuleID());
 
          for (typename container_type::const_iterator mod = zip.begin(); mod != zip.end(); ++mod) {
 
            const double Tmin  = root->getT() + mod->Tmin_ns;
            const double Tmax  = root->getT() + mod->Tmax_ns;
 
            JFrame_t&    frame = clone[mod->first];
 
            if (!frame.empty()) {
 
              frame.fast_forward(root->getT() - TMaxEvent_ns);
 
              for (typename JFrame_t::const_iterator i = frame.get(); *i <= Tmax; ++i) {
 
                if (*i >= Tmin) {
                  buffer.push_back(JHitR1(frame.getModuleID(),
                                          frame.getPosition(),
                                          frame.getJHit(*i)));
                }
              }
            }
          }
 
          if (distance          (buffer.begin(), buffer.end()) >= parameters.numberOfHits    - 1 &&
              getNumberOfModules(buffer.begin(), buffer.end()) >= parameters.numberOfModules - 1) {
 
            if (distance(buffer.begin(), buffer.end()) < parameters.factoryLimit) {
 
              JTriggerInput::iterator q = clusterize(buffer.begin(), buffer.end(), match, parameters.numberOfHits - 1);
 
              if (distance          (buffer.begin(), q) >= parameters.numberOfHits    - 1 &&
                  getNumberOfModules(buffer.begin(), q) >= parameters.numberOfModules - 1) {
 
                if (check(*root, buffer.begin(), q)) {
                  *out = JEvent(inputL1.getDAQChronometer(), *root, buffer.begin(), q, this->getTriggerBit());
                  ++out;
                }
              }
 
            } else {
 
              // Anomalous large event
        
              JEvent event(inputL1.getDAQChronometer(), *root, buffer.begin(), buffer.end(), this->getTriggerBit());
 
              event.addTriggerBit(FACTORY_LIMIT);
 
              *out = event;
              ++out;
            }
          }
        }
      }
    }  
 
 
    /**
     * Check event.
     *
     * \param  root             root hit
     * \param  __begin          begin of data
     * \param  __end            end   of data
     * \return                  true if accepted; else false
     */
    template<class T>
    inline bool check(const JHitR1& root, T __begin, T __end) const
    {
      using namespace std;
      using namespace JPP;
 
      if (!this->empty()) {
 
        const double U = getIndexOfRefraction() * getInverseSpeedOfLight();
 
        for (const_iterator u = this->begin(); u != this->end(); ++u) {
 
          const JVector3D p0 = root.getPosition() + *u;               // vertex position
          const double    t0 = root.getT()  -  u->getLength() * U;    // vertex time
 
          size_t   n = 1;                                             // number of hits
          set<int> m({root.getModuleID()});                           // number of modules
 
          for (T p = __begin; p != __end; ++p) {
 
            const double t1 = p->getT()  -  p0.getDistance(*p) * U;   // vertex time
 
            if (fabs(t1 - t0) <= parameters.TMaxExtra_ns) {
              n += 1;
              m.insert(p->getModuleID());
            }
          }
 
          if (n        >= (size_t) parameters.numberOfHits    &&
              m.size() >= (size_t) parameters.numberOfModules) {
            return true;
          }
        }
 
        return false;
 
      } else {
 
        return true;
      }
    }
 
  private:
    JParametersHelper<JParameters> parameters;
    JModuleMapper_t                mapper;
    mutable std::vector<JHitR1>    buffer;
    JMatch3G<JHitR1>               match;
    double                         TMaxEvent_ns;
  };
}
 
#endif