SlowMMTrigger_module.cc

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////
// Class:       SlowMMTrigger
// Module Type: filter
// File:        SlowMMTrigger_module.cc
//
// Generated at Mon Aug 26 13:26:30 2013 by Zukai Wang using artmod
// from cetpkgsupport v1_02_00.
////////////////////////////////////////////////////////////////////////

//ART Framework Include
#include "art/Framework/Core/EDFilter.h"
#include "art/Framework/Core/ModuleMacros.h"
#include "art/Framework/Principal/Event.h"
#include "art/Framework/Principal/Handle.h"
#include "art/Framework/Principal/Run.h"
#include "fhiclcpp/ParameterSet.h"
#include "art/Framework/Services/Optional/TFileService.h"       
#include "art/Framework/Services/Optional/TFileDirectory.h"
#include "canvas/Persistency/Common/FindOneP.h"
#include "canvas/Persistency/Common/Assns.h"
#include "canvas/Persistency/Common/FindManyP.h"

// NOvADDT includes
#include "DDTBaseDataProducts/BaseProducts.h"
#include "DDTBaseDataProducts/Boundary.h"
#include "DDTBaseDataProducts/BoundaryList.h"
#include "DDTBaseDataProducts/CompareDAQHit.h"
#include "DDTBaseDataProducts/HitList.h"
#include "DDTBaseDataProducts/DAQHit.h"
#include "DDTBaseDataProducts/GroupedHitList.h" 
#include "DDTBaseDataProducts/TriggerDecision.h"
#include "DDTBaseDataProducts/Track.h"
#include "DDTBaseDataProducts/Track3D.h"
#include "DDTUtilities/AssociationUtil.h"

#include <algorithm> 
#include <fstream>
#include <time.h>
#include <fstream>

namespace novaddt {
  class SlowMMTrigger;
}

class novaddt::SlowMMTrigger : public art::EDFilter {
public:
  explicit SlowMMTrigger(fhicl::ParameterSet const & p);
  virtual ~SlowMMTrigger();

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

  //This function should return an int showidng the district of the hit: 
  //0 -> non surface
  bool SurfAssign(DAQHit const & hptr);
  
  //Return the indices of the earliest and latest matching elements in the other view:
  bool ViewCheck(unsigned p1, TDC::value_type T1,
                 unsigned p2, TDC::value_type T2); 

  //This function return the index of the hit which closest to th time you assigned
  unsigned IndexLocator(HitList & hitlist,
                        TDC::value_type T);

  
  float Dist_(float c1, float c2,
              float p1, float p2);

private:
  unsigned _prescale;

  std::string _hitslabel;
  std::string _instance;

  float _hf_min;     //The threshold of hits/distance

  unsigned _deltaT;    //Time resolution (TDC) for single hit
  unsigned _sigmaT;
  unsigned _minHits;   //precut on hits number limit
  float _deltP;    //This is used for surface determination and also for view matching determination
  float _deltC;    //For surface determination
  float _sigmaL;    // Triangle cut for the max additional ratio allowed
  float _Pmin;
  float _Pmax;
  float _Cmax;
  float _DistMin;   //Min distance requirement    
  float _vproMin;   //Max beta: 0.01
  float _vproMax;   //Min beta: 4e-5
  unsigned _trigger_counts;
  float _df;        //CCB mechanism parameter 
  int   _CPUticks;  //If you consume more than this amount of CPU time, just trigger the lowest safe timewindow
};

novaddt::SlowMMTrigger::SlowMMTrigger(fhicl::ParameterSet const & p)
  :
  _prescale   (p.get< unsigned    >("prescale")),
  _hitslabel  (p.get< std::string >("hits_label")),
  _instance   (p.get< std::string >("instance"  )),
  _hf_min     (p.get< float    >("hf_min")),     //Min hit number density in unit distance
  _deltaT     (128),
  _sigmaT     (64),                         //Time difference allowence when attaching hits online & index locator
  _minHits    (p.get<unsigned>("minHits")),
  _deltP      (p.get<float>("deltP")),
  _deltC      (p.get<float>("deltC")),
  _sigmaL     (p.get<float>("sigmaL")),
  _Pmin       (p.get<float>("Pmin")), 
  _Pmax       (p.get<float>("Pmax")), 
  _Cmax       (p.get<float>("Cmax")),            
  _DistMin    (p.get<float>("DistMin")),
  _vproMin    (p.get<float>("vproMin")),
  _vproMax    (p.get<float>("vproMax")),
  _trigger_counts(0),
  _df         (0.1),
  _CPUticks   (p.get<int>("CPUticks"))
{
  // Call appropriate Produces<>() functions here.
  produces<std::vector<TriggerDecision>>();
}

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

bool novaddt::SlowMMTrigger::filter(art::Event & e)
{
  //Record the starting time:
  clock_t start_T = clock();
  bool result = false;
  std::unique_ptr<std::vector<TriggerDecision>> 
    trigger_decisions(new std::vector<TriggerDecision>);

  //Get the hits:
  art::Handle<novaddt::HitList> hits;
  e.getByLabel(_hitslabel, _instance, hits);

  HitList hl[2];
  //Seperate the hits by view, preserve the time order:
  for (unsigned i = 0; i!= hits->size(); ++i ) {
    unsigned v = 0;
    if ( ((*hits)[i]).View().val==daqchannelmap::Y_VIEW ) v = 1;
    hl[v].emplace_back((*hits)[i]);
  }

  //  std::cout<<"Xhits: "<<hl[0].size()<<"   Yhits: "<<hl[1].size();
  if (hl[0].size()<_minHits+2 || hl[1].size()<_minHits+2) {
    e.put(std::move(trigger_decisions));
    return false;
  }
  
  //Only record X hits
  HitList c_PS;     //They will contain all Valid non surface hits on each view
  HitList s_PS;     //They will contain all Valid surface and mate hits on each view

  //View check
  //Goal: doing view check for each X hit, and preserve the time order
  for (unsigned i = 0; i!= hl[0].size(); ++i ) {
    DAQHit const& hx = hl[0][i];
    TDC::value_type Tx = hx.TDC().val;
    unsigned p1= hx.Plane().val;
    //Do index locator calculation to locate the indices range:
    unsigned start_id = IndexLocator(hl[1], Tx-_deltaT);
    unsigned stop_id  = IndexLocator(hl[1], Tx+_deltaT); 
    bool met = false;
    //Look for your mate, once found, push back and break:
    for (unsigned j = start_id; j < stop_id; ++j) {
      DAQHit const& hy = hl[1][j];
      unsigned p2 = hy.Plane().val;
      TDC::value_type Ty = hy.TDC().val;
 
      if (ViewCheck(p1,Tx,p2,Ty) ) {
        if (SurfAssign(hx) || SurfAssign(hy)) {
          s_PS.emplace_back(hx);
          break;
        }
        else if (!met) c_PS.emplace_back(hx);
        met = true;
      }
    }
  }
  /*
  std::cout<<"  ViewCheckTicks: "<<clock()-start_T;
  std::cout<<"  SurfXhits: "<<s_PS.size()
           <<"  ContainedXhits: "<<c_PS.size();
  */
  
  unsigned x_con = c_PS.size();
  if (x_con < _minHits+2 || s_PS.size()<2) {
    e.put(std::move(trigger_decisions));
    return false;
  }

  //Making possible combinations of entry & exit hits:
  //Hints: you usually get more y hits than x hits, so look for candidates in X view 
  for (unsigned i=0; i<s_PS.size()-1 && !result; ++i) {
    DAQHit const& entry = s_PS[i]; 
    TDC::value_type  Ti = entry.TDC().val;
    unsigned pi = entry.Plane().val;
    unsigned ci = entry.Cell().val; 
    unsigned ini = IndexLocator(c_PS, Ti)+1;

    if (ini+_minHits+1 > x_con ) break;  //Don't go to the next hit, you won't get enough hits

    //Notice: the entry hit and the exit hit may not be in the same view:
    //Notice: exit hit should be looping from back to prevent double triggering!
    //According to the chronological order of hits, once j is triggered, 
    //you don't have to go down further, later triggers will be included in the time window!
    for (unsigned j=s_PS.size()-1; j>i && !result; --j) {
      DAQHit const&  exit = s_PS[j]; 
      TDC::value_type Tf  = exit.TDC().val;
      
      unsigned pf = exit.Plane().val;
      unsigned cf = exit.Cell().val;
      
      if (Tf+_deltaT < c_PS[ini+1+_minHits].TDC().val )  break;
      else if (Tf > c_PS[x_con-1].TDC().val+_deltaT )  continue;
      
      //Set the plane range:
      unsigned p_max = std::max(pi,pf);
      unsigned p_min = std::min(pi,pf);
      //Set the cell range:
      unsigned c_max = std::max(ci,cf);
      unsigned c_min = std::min(ci,cf);

      //If you are not traveling through enough planes, you also get cut:
      if (p_max < p_min+2*_deltP) continue;

      float Distance = Dist_(ci,cf,pi,pf);
      float vpro = (Tf-Ti)/Distance; // Always positive since hits are sorted
 
      if (Distance < _DistMin) continue;
      if (Tf > Ti + unsigned(_vproMax*Distance)+_sigmaT ) continue;
      if (Tf+_sigmaT < Ti + unsigned(_vproMin*Distance) ) continue;
      
      unsigned fin =  IndexLocator(c_PS, Tf);
      
      if (fin<ini+Distance*_hf_min) continue;

      //Before looping the mid point one by one, check if you are running out of time,
      if (clock()-start_T > _CPUticks) {
        result = true;
        trigger_decisions->
        emplace_back(Ti, Tf - Ti, 
                       daqdataformats::trigID::TRIG_ID_DATA_SLOWMONO, _prescale);

        /*      
        std::cout<<"  ot: "<<1<<"  ticks: "<<clock()-start_T
          <<"  DT: "<<Tf - Ti<<std::endl;
        */      
        e.put(std::move(trigger_decisions));
        return result;
      } 

      
      //Hmm, here comes the meat of this piece of code:
      //Loop through the non-boundary hits in the view, in the given range:
      unsigned h_counts = 0;
 
      float disCut   = Distance*_sigmaL;
      unsigned h_cut = Distance*_hf_min;        

      //Here we implement the CCB(continuity-check and break) mechanisim
      //First, calculate the time window for checking the previous hit on track:
      unsigned _dP  = (p_max-p_min)*_df;
      if (_dP<2) _dP=2;
      else if (_dP>4) _dP=4;
      float dt = (Tf-Ti)*_df;
      if (dt>256) dt=256;  //4us, enough for beta-4 to travel 12 cm!
      TDC::value_type T_last  = Ti;
      unsigned P_last = pi;      
      unsigned C_last = ci;

      for (unsigned k=ini; k!=fin;++k) {
        DAQHit const& midP = c_PS[k];
        unsigned pk = midP.Plane().val;
        unsigned ck = midP.Cell().val; 
        TDC::value_type Tk  = midP.TDC().val;
        if (pk>p_max) continue;
        if (pk<p_min) continue;
        if (ck>c_max) continue;
        if (ck<c_min) continue;

        //Now you can check whether it is going forward:
        if (pf>pi) {
          if ( pk> P_last+_dP || pk < P_last){
            if (Tk-T_last>dt) break;
            else continue;
          }

          if (cf > ci && ck < C_last ) continue;
          if (cf < ci && ck > C_last ) continue;
        }

        else {
          if (pk+_dP < P_last || pk > P_last ) {
            if (Tk-T_last>dt) break;
            else continue;
          }
          if (pk+_dP < P_last || pk > P_last || ck > C_last+_deltC || ck+_deltC<C_last) continue;
          if (cf > ci && ck < C_last ) continue;
          if (cf < ci && ck > C_last ) continue;
        }
                
        float d1 = Dist_(ci,ck,pi,pk);
        float d2 = Dist_(cf,ck,pf,pk);

        if (d1+d2 > Distance+disCut) continue;
        float t1 = Tk-Ti;   // From start to mid
        float t2 = Tf-Tk;   // From mid to end

        if (d1>d2) {
          if (t1 < d1*vpro + _sigmaT && t1 + _sigmaT > d1*vpro) {
            ++h_counts;
          //refresh the previous hit:
            P_last = pk;
            C_last = ck;
            T_last = Tk;
          }
        }

        else {
          if (t2 < d2*vpro + _sigmaT && t2 + _sigmaT > d2*vpro ) {
            ++h_counts;
            //refresh the previous hit:
            P_last = pk;
            C_last = ck;
            T_last = Tk;
          }
        }                 
      }

      if (h_counts>h_cut && h_counts>_minHits-2) {
        result = true;  
        trigger_decisions->
          emplace_back(Ti-2*_deltaT, Tf - Ti+4*_deltaT, 
                         daqdataformats::trigID::TRIG_ID_DATA_SLOWMONO, _prescale);
        /*
        std::cout<<"  ot: "<<0<<"  ticks: "<<clock()-start_T
                 <<"  DT: "<<Tf - Ti+4*_deltaT;
        */
        break;
      }
    }
  }     
  /*
  std::cout<<"  ot: "<<0<<"  ticks: "<<clock()-start_T
           <<"  DT: "<<-1<<std::endl;
  */
  e.put(std::move(trigger_decisions));
  return result;

}

bool novaddt::SlowMMTrigger::SurfAssign(DAQHit const & hit) 
{
  //Check if this is close to surface
  if      ( hit.Cell().val  > _Cmax - _deltC   ) return true;
  else if ( hit.Cell().val  < _deltC           ) return true;
  else if ( hit.Plane().val < _deltP           ) return true;
  else if ( hit.Plane().val > _Pmax - _deltP   ) return true;

  return false;
}

//Note: Ts < Tm (Looking forwards). Time check is already done in the main loop.
bool novaddt::SlowMMTrigger::ViewCheck(unsigned p1, TDC::value_type T1,
                                       unsigned p2, TDC::value_type T2)
{
  if (T1>T2+_deltaT) return false;
  if (T2>T1+_deltaT) return false;
  if (p1 > p2 ) {
    if (p1-p2 < _deltP) return true;
  }
  else if (p2-p1 < _deltP) return true;

  return false; 
}

//It should return the index before T
unsigned novaddt::SlowMMTrigger::IndexLocator(HitList & hl,
                                              TDC::value_type T)
{
  unsigned i = 0;
  unsigned f = hl.size()-1;

  if (T+_deltaT < hl[i].TDC().val ) return i; 
  if (T> hl[f].TDC().val+_deltaT  ) return f;

  while (f>i+1) {
    unsigned m = (f+i)/2;
    if (T >= hl[m].TDC().val ) i = m;
    else f = m; 
  }
  return i;
}


float novaddt::SlowMMTrigger::Dist_(float c1, float c2,
                                    float p1, float p2)
{ 
  return sqrt((p1-p2)*(p1-p2)*2.82072+(c1-c2)*(c1-c2));
}


DEFINE_ART_MODULE(novaddt::SlowMMTrigger)