NNFastMMTrigger_module.cc

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////
// Class:       NNFastMMTrigger
// Module Type: filter
// File:        NNFastMMTrigger_module.cc
//
// Created at Nov 22 2015 by Enhao Song 
////////////////////////////////////////////////////////////////////////

#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 "art/Framework/Principal/Run.h"
#include "art/Framework/Principal/SubRun.h"
#include "art/Framework/Services/Optional/TFileService.h"
#include "canvas/Utilities/InputTag.h"
#include "fhiclcpp/ParameterSet.h"
#include "fhiclcpp/ParameterSetRegistry.h"
#include "messagefacility/MessageLogger/MessageLogger.h"

#include "DDTUtilities/AssociationUtil.h"
#include "DDTBaseDataProducts/BaseProducts.h"
#include "DDTBaseDataProducts/Boundary.h"
#include "DDTBaseDataProducts/BoundaryList.h"
#include "DDTBaseDataProducts/HitList.h"
#include "DDTBaseDataProducts/DAQHit.h"
#include "DDTBaseDataProducts/GroupedHitList.h"
#include "DDTBaseDataProducts/TriggerDecision.h"
#include "DDTBaseDataProducts/EventHeader.h"
#include "DDTBaseDataProducts/Track.h"  

#include <string>
#include <vector>
#include <algorithm>
#include <numeric>
#include <cmath>
#include <map>

#include "ExoticsTriggers/KerasModel.h"

namespace novaddt {
  class NNFastMMTrigger;
}
class novaddt::NNFastMMTrigger : public art::EDFilter {
  public:
    explicit NNFastMMTrigger(fhicl::ParameterSet const & p);
    virtual ~NNFastMMTrigger();
    virtual bool filter(art::Event & e) override;
    void beginJob() override;
    void endJob() override;
    void StdevCellsPerPlane(HitList const &hits, double &stdev_cells_per_plane_xview, double &stdev_cells_per_plane_yview ) const;
    void NumberOfCellsPerLength(HitList const &hits,double tracklength_of_xview, double tracklength_of_yview, double &number_of_cells_per_length_xview,double &number_of_cells_per_length_yview) const;
    float Distance(unsigned const &xcell1, unsigned const  &ycell1, unsigned const &plane1,unsigned const &xcell2, unsigned const &ycell2 , unsigned const &plane2) const;
    bool SurfAssign(DAQHit const & hit) const;
    double WeightedCenterCut(HitList const &hits,float PX_min,float PX_max,float PY_min,float PY_max,float CX_min, float CX_max, float CY_min, float CY_max, double &weighted_off_center_xx,double &weighted_off_center_xz,double &weighted_off_center_yy,double &weighted_off_center_yz) const;
    double NumberOfHitsInOverlapPlanesCut(HitList const &hits,float PX_min,float PX_max,float PY_min,float PY_max ) const;
    int NumberOfSurfaceHits(HitList const &hits) const;
  private:
    unsigned _prescale;
    std::string _slicelabel,_sliceinstance, _monopole_id_lib_path;
    int _xMin, _xMax, _yMin, _yMax, _zMin, _zMax, _deltx, _delty, _deltz, _sigmaP;
    int PX_max, PY_max, PX_min, PY_min, CX_max, CX_min, CY_max, CY_min;
    TDC::value_type TX_max, TY_max, TX_min, TY_min;
    int _verbose; 
    unsigned _trigger_counts; 
    keras::KerasModel * _model=NULL;
    void LinFit(const std::vector<double>& x, const std::vector<double>& y, double *fitpar);
};
novaddt::NNFastMMTrigger::NNFastMMTrigger(fhicl::ParameterSet const & p)
  :
    _prescale(p.get<unsigned>("prescale")),
    _slicelabel(p.get<std::string>("slice_label")),
    _sliceinstance(p.get<std::string>("slice_instance")),
    _monopole_id_lib_path(p.get<std::string>("MONOPOLEIDLibPath")),
    _xMin(p.get<int>("x_min")),
    _xMax(p.get<int>("x_max")),
    _yMin(p.get<int>("y_min")),
    _yMax(p.get<int>("y_max")),
    _zMin(p.get<int>("z_min")),
    _zMax(p.get<int>("z_max")),
    _deltx(p.get<int>("x_delt")),
    _delty(p.get<int>("y_delt")),
    _deltz(p.get<int>("z_delt")),
    _sigmaP(p.get<int>("sigma_P")),
    _verbose(p.get<int>("verbose")),
    _trigger_counts(0)
{
  // Call appropriate Produces<>() functions here.
  produces<std::vector<TriggerDecision>>();
}

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

void novaddt::NNFastMMTrigger::beginJob(){
    _model = new keras::KerasModel(_monopole_id_lib_path);
}
void novaddt::NNFastMMTrigger::endJob(){
  delete _model;
}

bool novaddt::NNFastMMTrigger::filter(art::Event & e)   {
  bool result = false;
  std::unique_ptr<std::vector<TriggerDecision>> 
    trigger_decisions(new std::vector<TriggerDecision>);
  art::Handle<std::vector<novaddt::HitList> > slices;
  e.getByLabel(_slicelabel, _sliceinstance, slices);

  for (unsigned i = 0; i!= slices->size(); ++i) {
    HitList const &hits = (*slices)[i];
    PX_max = 0;//P is short for plane.
    CX_max = 0;//C is short for cell. 
    CY_max = 0;
    PY_max = 0;
    PX_min = 99999;
    CX_min = 99999;
    CY_min = 99999;
    PY_min = 99999;
    TX_max = hits.front().TDC().val;
    TY_max = TX_max;
    TX_min = hits.back().TDC().val;
    TY_min = TX_min;
    if(_verbose != 0){ 
      int nb_hits = hits.size();
      std::cout<<"slice number = "<<i+1<<" number of hits: "<<nb_hits<<std::endl;
    }
    int  entry_id = 0, nx = 0, ny = 0;
    bool PassThrough = false, Intrusion   = false;
    double meanADC = 0, stdevADC = 0, sq_sum_ADC = 0;

    std::vector<double> x_hit, y_hit, zy_hit, zx_hit;
    for (auto const hit: hits) {
      if(_verbose != 0) std::cout<<"TDC: "<<hit.TDC().val << "    ADC: "<<hit.ADC().val<<"   Plane: "<<hit.Plane().val<<"   Cell:  "<<hit.Cell().val<<std::endl;
      //X view hits
      meanADC += hit.ADC().val;
      sq_sum_ADC += hit.ADC().val*hit.ADC().val;
      if (hit.View().val == daqchannelmap::X_VIEW) {
        if(_verbose != 0) std::cout<<"X hits"<<std::endl;
        ++nx;
        x_hit.push_back(hit.Cell().val);
        zx_hit.push_back(hit.Plane().val);
        if (PX_max < hit.Plane().val)   PX_max = hit.Plane().val;
        if (PX_min > hit.Plane().val)   PX_min = hit.Plane().val;
        if (CX_max < hit.Cell().val)    CX_max = hit.Cell().val;
        if (CX_min > hit.Cell().val)    CX_min = hit.Cell().val;
        if (TX_max < hit.TDC().val)     TX_max = hit.TDC().val;
        if (TX_min > hit.TDC().val)     TX_min = hit.TDC().val;  
        if (!Intrusion) {
          if (hit.Cell().val < _xMin+_deltx) {
            entry_id  = 1;
            Intrusion = true;
          }
          else if (hit.Cell().val> _xMax-_deltx) {
            entry_id  = 2;
            Intrusion = true;
          }
        }
        else if (!PassThrough) {
          if (hit.Cell().val < _xMin+_deltx && entry_id != 1) PassThrough = true;
          else if (hit.Cell().val> _xMax-_deltx && entry_id != 2) PassThrough = true;
        }
      } 
      //Y view hits
      else {
        if(_verbose != 0)  std::cout<<"Y hits"<<std::endl;
        ++ny;
        y_hit.push_back(hit.Cell().val);
        zy_hit.push_back(hit.Plane().val);
        if (PY_max < hit.Plane().val)   PY_max   = hit.Plane().val;
        if (PY_min > hit.Plane().val)   PY_min   = hit.Plane().val;
        if (CY_max < hit.Cell().val)    CY_max  = hit.Cell().val;
        if (CY_min > hit.Cell().val)    CY_min  = hit.Cell().val;
        if (TY_max < hit.TDC().val )   TY_max   = hit.TDC().val;
        if (TY_min > hit.TDC().val )   TY_min   = hit.TDC().val;
        if (!Intrusion) {
          if (hit.Cell().val < _yMin+_delty) {
            entry_id  = 3;
            Intrusion = true;
          }
          else if (hit.Cell().val> _yMax-_delty) {
            entry_id  = 4;
            Intrusion = true;
          }
        }
        else if (!PassThrough) {
          if (hit.Cell().val < _yMin+_delty && entry_id != 3) PassThrough = true;
          else if (hit.Cell().val> _yMax-_delty && entry_id != 4) PassThrough = true;
        }
      }
      //No matter which view:
      if (!Intrusion) {
        if (hit.Plane().val < _zMin+_deltz) {
          entry_id  = 5;
          Intrusion = true;
        }
        else if (hit.Plane().val > _zMax-_deltz) {
          entry_id  = 6;
          Intrusion = true;
        }
      }
      else if (!PassThrough) {
        if (hit.Plane().val < _zMin+_deltz && entry_id != 5) PassThrough=true;
        else if (hit.Plane().val > _zMax-_deltz && entry_id != 6) PassThrough=true;
      }
    }
    meanADC = meanADC/(nx+ny);

    //Min-hits on both view:
    if(_verbose != 0)       std::cout<<"xhits number : "<<nx<<"   yhits number : "<<ny<<std::endl;
    if (nx<10 || ny<10) continue;
    if(_verbose != 0)       std::cout<<"minimum hits requirement met!"<<std::endl;
    //Passing Through Test:
    if (!PassThrough) continue;
    if(_verbose != 0)       std::cout<<"through!"<<std::endl;
    if(_verbose != 0)       std::cout<<"X view Plane range: "<<PX_min<<" "<<PX_max<<std::endl;
    if(_verbose != 0)       std::cout<<"X view Cell range: "<<CX_min<<" "<<CX_max<<std::endl;
    if(_verbose != 0)       std::cout<<"X view time range: "<<TX_min<<" "<<TX_max<<std::endl;

    if(_verbose != 0)       std::cout<<"Y view Plane range: "<<PY_min<<" "<<PY_max<<std::endl;
    if(_verbose != 0)       std::cout<<"Y view Cell range: "<<CY_min<<" "<<CY_max<<std::endl;
    if(_verbose != 0)       std::cout<<"Y view time range: "<<TY_min<<" "<<TY_max<<std::endl;

    //Matching Test: 
    double deltaP = std::max(PX_max,PY_max) - std::min(PX_min,PY_min);
    double deltaTDC = std::max(TX_max,TY_max) - std::min(TX_min,TY_min);
    if (std::min(PX_max,PY_max) < std::max(PX_min,PY_min)+_sigmaP) continue; //Require minimum overlapping planes
    if (std::min(TX_max,TY_max) < std::max(TX_min,TY_min) ) continue;
    //double preselection = 1;

    unsigned dT = std::max(TX_max,TY_max)-std::min(TX_min,TY_min);

    float tracklength =Distance(CX_max, CY_max, PX_max, CX_min, CY_min, PX_min); 
    double number_of_surface_hits = NumberOfSurfaceHits(hits);
    if(_verbose != 0)       std::cout<<"number of surface hits: "<<number_of_surface_hits<<std::endl;
    if(_verbose != 0)       std::cout<<"slice diagonal length: "<<tracklength<<std::endl;
    if(_verbose != 0)       std::cout<<"number of hits divided by slice diagonal length"<<hits.size()/tracklength<<std::endl;

    // Linear fit of the hits in XZ/YZ views.
    double fitpar[3];
    LinFit(x_hit, zx_hit, fitpar);
    double r2x = 0;
    if( std::isfinite(fitpar[2]))
      r2x = fitpar[2];
    LinFit(y_hit, zy_hit, fitpar);
    //double r2y = 0;

    double number_of_cells_per_length_xview=0;
    double number_of_cells_per_length_yview=0;
    //double tracklength_of_xview=Distance(CX_max, 0, PX_max, CX_min, 0, PX_min);
    //double tracklength_of_yview=Distance(0, CY_max, PY_max, 0, CY_min, PY_min);
//    NumberOfCellsPerLength(hits,tracklength_of_xview,tracklength_of_yview,number_of_cells_per_length_xview,number_of_cells_per_length_yview);
    if(_verbose != 0) std::cout<<"number of cells per length xview: "<<number_of_cells_per_length_xview<<" yview: "<<number_of_cells_per_length_yview<<std::endl;
    double stdev_cells_per_plane_xview=0, stdev_cells_per_plane_yview=0;
    StdevCellsPerPlane(hits, stdev_cells_per_plane_xview, stdev_cells_per_plane_yview);
    if(_verbose != 0) std::cout<<"stdev cells per plane xview: "<<stdev_cells_per_plane_xview<<" yview: "<<stdev_cells_per_plane_yview<<std::endl;
    double weighted_off_center_xx=0, weighted_off_center_xz=0, weighted_off_center_yy=0, weighted_off_center_yz=0;
    WeightedCenterCut(hits, PX_min, PX_max, PY_min, PY_max, CX_min,  CX_max,  CY_min,  CY_max, weighted_off_center_xx, weighted_off_center_xz, weighted_off_center_yy, weighted_off_center_yz );
    if(_verbose != 0) std::cout<<"weighted off center xz: "<<weighted_off_center_xz<<" yz: "<<weighted_off_center_yz<<std::endl;
    //double percentage_of_hits_in_overlap_planes = NumberOfHitsInOverlapPlanesCut(hits,PX_min,PX_max,PY_min,PY_max);
    //if(_verbose != 0) std::cout<<"percent of hits in overlap planes: "<<percentage_of_hits_in_overlap_planes<<std::endl;
    stdevADC = std::sqrt( sq_sum_ADC*1.0/(nx+ny) - meanADC*meanADC);
    //Keras model
    keras::DataChunk *sample = new keras::DataChunkFlat();
    //Init input to keras model with normalization.
    float myfloats[] = {(hits.size()/tracklength-0.0206254683435)/94.6383405961,
    (stdev_cells_per_plane_xview-0.816496580928)/12.9743518339,
    (stdev_cells_per_plane_yview-0.816496580928)/15.2366123639,
    (weighted_off_center_yy-0.0)/0.499630220584,
    (number_of_surface_hits-1.0)/30304.0,
    (meanADC-239.074074074)/3603.61013645,
    (stdevADC-8.75437337874)/1912.00668817,
    (ny-3.0)/78720.0,
    (deltaTDC-0.0)/27478.0,
    (deltaP-3.0)/892.0,
    (r2x-0.0)/1.0 };

    std::vector<float> v (myfloats, myfloats + sizeof(myfloats) / sizeof(float) );
    if(_verbose!=0){
      for (auto data_i : v)
        std::cout<<data_i<<std::endl;
    }
    sample->set_data(v);
    double score = _model->compute_output(sample)[0];
    if(_verbose != 0)  std::cout<< "Score is: "<< score<< std::endl;
    delete sample;
    if(score<0.5) continue;
    ++_trigger_counts;

    if (_trigger_counts%_prescale == _prescale-1) {
      trigger_decisions->
        emplace_back(std::min(TX_min,TY_min), 
            dT,  
            daqdataformats::trigID::TRIG_ID_DATA_FASTMONONN, _prescale);
      result = true;
    }
  }
  e.put(std::move(trigger_decisions));
  return result;
}

void novaddt::NNFastMMTrigger::LinFit(const std::vector<double>& x_val, const std::vector<double>& y_val, double *fitpar){
  // http://en.wikipedia.org/wiki/Simple_linear_regression
  const auto n  = x_val.size();
  const auto x  = (std::accumulate(x_val.begin(), x_val.end(), 0.0))/n;                       // <x>
  const auto y  = (std::accumulate(y_val.begin(), y_val.end(), 0.0))/n;                       // <y>
  const auto xx = (std::inner_product(x_val.begin(), x_val.end(), x_val.begin(), 0.0))/n;     // <xx>
  const auto yy = (std::inner_product(y_val.begin(), y_val.end(), y_val.begin(), 0.0))/n;     // <yy>
  const auto xy = (std::inner_product(x_val.begin(), x_val.end(), y_val.begin(), 0.0))/n;     // <xy>

  const auto b  = (xy - x*y)/(xx - x*x);                                                    // slope
  const auto a  = y - b*x;                                                                  // intercept
  const auto r  = (xy - x*y)/sqrt((xx - x*x)*(yy - y*y));                                   // Rxy - coeffcient of determination
  fitpar[0] = a;
  fitpar[1] = b;
  fitpar[2] = r*r;
}

void novaddt::NNFastMMTrigger::NumberOfCellsPerLength(HitList const &hits,double tracklength_of_xview, double tracklength_of_yview, double &number_of_cells_per_length_xview,double &number_of_cells_per_length_yview) const{
  std::vector<std::vector<int>> x_hits_plane(PX_max-PX_min+1,std::vector<int>(0));
  std::vector<std::vector<int>> y_hits_plane(PY_max-PY_min+1,std::vector<int>(0)); 
  for (auto const hit: hits) {
    if (hit.View().val == daqchannelmap::X_VIEW) {
      //std::cout<<"X hits"<<std::endl;
      x_hits_plane[hit.Plane().val-PX_min].push_back(hit.Cell().val); 
    }
    if (hit.View().val == daqchannelmap::Y_VIEW) {
      //std::cout<<"Y hits"<<std::endl;
      y_hits_plane[hit.Plane().val-PY_min].push_back(hit.Cell().val); 
    }
  }
  //following is to get unique number of cells hit per plane
  std::vector<int> x_cells_plane(PX_max-PX_min+1,0);
  std::vector<int> y_cells_plane(PY_max-PY_min+1,0); 
  for (unsigned i=0;i!=x_hits_plane.size();i++){
    std::sort(x_hits_plane[i].begin(),x_hits_plane[i].end());
    std::vector<int>::iterator it;
    it = std::unique (x_hits_plane[i].begin(),x_hits_plane[i].end());
    x_hits_plane[i].resize(std::distance(x_hits_plane[i].begin(),it));
    x_cells_plane[i]=x_hits_plane[i].size();
  }
  for (unsigned i=0;i!=y_hits_plane.size();i++){
    std::sort(y_hits_plane[i].begin(),y_hits_plane[i].end());
    std::vector<int>::iterator it;
    it = std::unique (y_hits_plane[i].begin(),y_hits_plane[i].end());
    y_hits_plane[i].resize(std::distance(y_hits_plane[i].begin(),it));
    y_cells_plane[i]=y_hits_plane[i].size();
  }
  //following is to get number of cells per track length in x view and y view
  int sum_of_cells=0;
  for (unsigned i=0;i!=x_cells_plane.size();i++){
    sum_of_cells+=x_cells_plane[i];
  }
  number_of_cells_per_length_xview=sum_of_cells/tracklength_of_xview;
  sum_of_cells=0;
  for (unsigned i=0;i!=y_cells_plane.size();i++){
    sum_of_cells+=y_cells_plane[i];
  }
  number_of_cells_per_length_yview=sum_of_cells/tracklength_of_yview;

}
void novaddt::NNFastMMTrigger::StdevCellsPerPlane(HitList const &hits, double &stdev_cells_per_plane_xview, double &stdev_cells_per_plane_yview ) const{
  //fluctuation of number of cells of each plane check by Enhao
  std::vector<std::vector<int>> x_hits_plane(PX_max-PX_min+1,std::vector<int>(0));
  std::vector<std::vector<int>> y_hits_plane(PY_max-PY_min+1,std::vector<int>(0)); 
  for (auto const hit: hits) {
    if (hit.View().val == daqchannelmap::X_VIEW) {
      //std::cout<<"X hits"<<std::endl;
      x_hits_plane[hit.Plane().val-PX_min].push_back(hit.Cell().val); 
    }
    if (hit.View().val == daqchannelmap::Y_VIEW) {
      //std::cout<<"Y hits"<<std::endl;
      y_hits_plane[hit.Plane().val-PY_min].push_back(hit.Cell().val); 
    }
  }
  double sum    = 0.;
  double mean   = 0.;
  double sq_sum = 0.;
  double stdev  = 0.;
  //double min_stdev = 0.;
  //following is to get unique number of cells hit per plane
  std::vector<int> x_cells_plane(PX_max-PX_min+1,0);
  std::vector<int> y_cells_plane(PY_max-PY_min+1,0); 
  for (unsigned i=0;i!=x_hits_plane.size();i++){
    std::sort(x_hits_plane[i].begin(),x_hits_plane[i].end());
    std::vector<int>::iterator it;
    it = std::unique (x_hits_plane[i].begin(),x_hits_plane[i].end());
    x_hits_plane[i].resize(std::distance(x_hits_plane[i].begin(),it));
    x_cells_plane[i]=x_hits_plane[i].size();
  }
  for (unsigned i=0;i!=y_hits_plane.size();i++){
    std::sort(y_hits_plane[i].begin(),y_hits_plane[i].end());
    std::vector<int>::iterator it;
    it = std::unique (y_hits_plane[i].begin(),y_hits_plane[i].end());
    y_hits_plane[i].resize(std::distance(y_hits_plane[i].begin(),it));
    y_cells_plane[i]=y_hits_plane[i].size();
  }
  sum = std::accumulate(x_cells_plane.begin(),x_cells_plane.end(),0.0);
  mean = sum / x_cells_plane.size();
  sq_sum = std::inner_product(x_cells_plane.begin(),x_cells_plane.end(),x_cells_plane.begin(),0.0);
  stdev = std::sqrt(sq_sum / x_cells_plane.size()-mean * mean);
  stdev_cells_per_plane_xview=stdev/mean;
  sum = std::accumulate(y_cells_plane.begin(),y_cells_plane.end(),0.0);
  mean = sum / y_cells_plane.size();
  sq_sum = std::inner_product(y_cells_plane.begin(),y_cells_plane.end(),y_cells_plane.begin(),0.0);
  stdev = std::sqrt(sq_sum / y_cells_plane.size()-mean * mean);
  stdev_cells_per_plane_yview=stdev/mean;
}
float novaddt::NNFastMMTrigger::Distance( unsigned const &xcell1, unsigned const  &ycell1, unsigned const &plane1,unsigned const &xcell2, unsigned const &ycell2 , unsigned const &plane2) const {
  return std::sqrt((plane1-plane2)*(plane1-plane2)*2.82072+(xcell1-xcell2)*(xcell1-xcell2)+(ycell1-ycell2)*(ycell1-ycell2));
}

bool novaddt::NNFastMMTrigger::SurfAssign(DAQHit const &hit) const{
  //Check if this is close to surface
  if      ( hit.Cell().val  > _xMax - _deltx   ) return true;
  else if ( hit.Cell().val  < _deltx           ) return true;
  else if ( hit.Plane().val < _deltz           ) return true;
  else if ( hit.Plane().val > _zMax - _deltz   ) return true;

  return false;
}
int novaddt::NNFastMMTrigger::NumberOfSurfaceHits(HitList const &hits) const{
  int number_of_surface_hits=0;
  for(unsigned i=0;i!=hits.size();++i){
    if(SurfAssign(hits[i])==true) number_of_surface_hits++;
  }
  return number_of_surface_hits;
}
double novaddt::NNFastMMTrigger::WeightedCenterCut(HitList const &hits,float PX_min,float PX_max,float PY_min,float PY_max,float CX_min, float CX_max, float CY_min, float CY_max, double &weighted_off_center_xx,double &weighted_off_center_xz,double &weighted_off_center_yy,double &weighted_off_center_yz ) const{
  float Weighted_PX_Center=0.5,Weighted_CX_Center=0.5,Weighted_PY_Center=0.5,Weighted_CY_Center=0.5,SumOfXHitsADC=0,SumOfYHitsADC=0;
  for (auto const hit: hits) {
    if (hit.View().val == daqchannelmap::X_VIEW) {
      SumOfXHitsADC+=hit.ADC().val;
      Weighted_PX_Center+=hit.ADC().val*(hit.Plane().val-PX_min);
      Weighted_CX_Center+=hit.ADC().val*(hit.Cell().val-CX_min);
    }
    if (hit.View().val == daqchannelmap::Y_VIEW) {
      SumOfYHitsADC+=hit.ADC().val;
      Weighted_PY_Center+=hit.ADC().val*(hit.Plane().val-PY_min);
      Weighted_CY_Center+=hit.ADC().val*(hit.Cell().val-CY_min);
    }
  }
  if(PX_max!=PX_min)  Weighted_PX_Center=Weighted_PX_Center/SumOfXHitsADC/(PX_max-PX_min);
  if(PY_max!=PY_min)  Weighted_PY_Center=Weighted_PY_Center/SumOfYHitsADC/(PY_max-PY_min);
  if(CX_max!=CX_min)  Weighted_CX_Center=Weighted_CX_Center/SumOfXHitsADC/(CX_max-CX_min);
  if(CY_max!=CY_min)  Weighted_CY_Center=Weighted_CY_Center/SumOfYHitsADC/(CY_max-CY_min);
  weighted_off_center_xx=std::abs(Weighted_CX_Center-0.5);
  weighted_off_center_xz=std::abs(Weighted_PX_Center-0.5);
  weighted_off_center_yy=std::abs(Weighted_CY_Center-0.5); 
  weighted_off_center_yz=std::abs(Weighted_PY_Center-0.5);
  return std::max(std::max(std::abs(Weighted_PX_Center-0.5),std::abs(Weighted_PY_Center-0.5)),std::max(std::abs(Weighted_CX_Center-0.5),std::abs(Weighted_CY_Center-0.5)));
}

double novaddt::NNFastMMTrigger::NumberOfHitsInOverlapPlanesCut(HitList const &hits,float PX_min,float PX_max,float PY_min,float PY_max ) const{
  int max=std::min(PX_max,PY_max);
  int min=std::max(PX_min,PY_min);
  double number_hits_in_overlap_planes=0;
  for (auto const hit: hits) {
    if(hit.Plane().val<=max && hit.Plane().val>=min)    number_hits_in_overlap_planes++;
  }
  return number_hits_in_overlap_planes/hits.size();
}

DEFINE_ART_MODULE(novaddt::NNFastMMTrigger)