FmmTriggerEvd_module.cc

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////
// Class:       FmmTriggerEvd
// Module Type: producer
// File:        FmmTriggerEvd_module.cc
//
// Generated at 21/01/2016 by Enhao Song
// this works like fast monopole trigger and filter then  produce the 
// slice that triggered
////////////////////////////////////////////////////////////////////////
#include "art/Framework/Core/EDProducer.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 "MCCheater/BackTracker.h"
//#include "Calibrator/Calibrator.h"
#include "Simulation/Particle.h"
#include "Simulation/ParticleNavigator.h"
#include "art/Framework/Principal/SubRun.h"
#include "canvas/Utilities/InputTag.h"
#include "Utilities/func/MathUtil.h"
#include "fhiclcpp/ParameterSet.h"
#include "art/Framework/Services/Optional/TFileService.h"

#include "Geometry/Geometry.h"
#include "GeometryObjects/CellGeo.h"
#include "GeometryObjects/Geo.h"
#include "GeometryObjects/PlaneGeo.h"

#include "NovaDAQConventions/DAQConventions.h"
#include "RecoBase/CellHit.h"
#include "RecoBase/RecoHit.h"
#include "RecoBase/Cluster.h"
#include "RecoBase/Track.h"

#include "TVector3.h"
#include "TTree.h"
#include "TProfile.h"

#include <algorithm> 
#include <utility>
#include <cmath>
#include <iostream>

namespace zcl {
  class FmmTriggerEvd;
}

class zcl::FmmTriggerEvd : public art::EDProducer {
  public:
    explicit FmmTriggerEvd(fhicl::ParameterSet const & pset);
    virtual ~FmmTriggerEvd();

    void produce(art::Event & evt);

  private:
    std::string fClusterInput;    ///< Input folder from cluster reco
    //Each entry is corresponding to 1 reconstructed slice

    int    evtID, subRunID, RunID;

    int    _Nslices;
    bool   _penetrated;
    int    _nxHits, _nyHits;
    int PX_max, PY_max, PX_min, PY_min, CX_min, CX_max, CY_min, CY_max;
    double _dpx, _dpy, _dcx, _dcy,_dp;
    double _rmsADC, _meanADC;
    double _cx,  _cy,  _cz;
    double _cxW, _cyW, _czW;

    //added by Enhao
    double _number_of_surface_hits;
    double _number_of_hits_per_length;
    double _number_of_cells_per_length_xview,_number_of_cells_per_length_yview;
    double _stdev_cells_per_plane_xview,_stdev_cells_per_plane_yview;
    double _weighted_off_center_xx,_weighted_off_center_xz,_weighted_off_center_yy,_weighted_off_center_yz;
    double _cxz,_cyz;
    double _cxzW,_cyzW;
    double _percent_of_hits_in_overlap_planes;


    double _tmin, _tmax;  //The time window of the slice

    const int _xMax       = 383; 
    const int _yMax       = 383; 
    const int _zMax       = 895; 
    const int _xdelt      = 35;
    const int _ydelt      = 35;
    const int _zdelt      = 15;

    double RMS(std::vector<double> & v, double & mean);
    bool IsSurfaceHit(int plane, int cell) const;
    void NumberOfCellsPerLength(rb::Cluster const &slice,double tracklength_of_xview, double tracklength_of_yview, double &number_of_cells_per_length_xview,double &number_of_cells_per_length_yview)const;
    void StdevCellsPerPlane(rb::Cluster const &slice, double &stdev_cells_per_plane_xview, double &stdev_cells_per_plane_yview ) const;
    double NumberOfHitsInOverlapPlanesCut(rb::Cluster const &slice,float PX_min,float PX_max,float PY_min,float PY_max ) const;
};

zcl::FmmTriggerEvd::FmmTriggerEvd(fhicl::ParameterSet const & pset) :
  fClusterInput      (pset.get< std::string > ("ClusterInput")      )
{
    produces< std::vector<rb::Cluster> >();
}

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

void zcl::FmmTriggerEvd::produce(art::Event & evt)
{
  evtID    = evt.id().event();
  subRunID = evt.subRun();
  RunID    = evt.run();

  art::ServiceHandle<geo::Geometry> geom;
  //  art::ServiceHandle<calib::Calibrator> calib;

  //=====================================================================================
  //Now, start pseudo trigger:

  double TX_max, TY_max, TX_min, TY_min;

  art::Handle< std::vector<rb::Cluster> > slices;
  evt.getByLabel(fClusterInput, slices);

    std::unique_ptr< std::vector<rb::Cluster> > clustercol(new std::vector<rb::Cluster>);

  _Nslices = slices->size(); 
  std::cout<<"number of slices: "<<_Nslices<<std::endl;

  //Instead of splitting the hits into 2 views, just directly take the entire 3D slice:
  for (unsigned sliceIdx = 0; sliceIdx != slices->size(); ++sliceIdx) {
    std::cout<<"slice ID: "<<sliceIdx<<std::endl;
    const rb::Cluster& slice = (*slices)[sliceIdx];
    std::vector<art::Ptr<rb::CellHit> > Xhitlist, Yhitlist;
    CX_min = 400;
    CX_max = 0;
    CY_min = 400;
    CY_max = 0;    
    PX_max = 0; 
    PY_max = 0;
    PX_min = 999;
    PY_min = 999;
    TX_max = -9e7;      
    TY_max = -9e7;
    TX_min = 9e9;
    TY_min = 9e9;

    _penetrated = false;
    _nxHits     = 0;
    _nyHits     = 0;
    _meanADC    = 0;
    _rmsADC     = 0;
    _cx         = 0;
    _cy         = 0;
    _cz         = 0;
    _cxW        = 0;
    _cyW        = 0;
    _czW        = 0;

    //added by enhao
    _number_of_surface_hits=0;
    _number_of_hits_per_length=0;
    _number_of_cells_per_length_xview=0;
    _number_of_cells_per_length_yview=0;
    _stdev_cells_per_plane_xview=0;
    _stdev_cells_per_plane_yview=0;
    _weighted_off_center_xx=0; 
    _weighted_off_center_xz=0; 
    _weighted_off_center_yy=0; 
    _weighted_off_center_yz=0; 
    _cxz=0;_cyz=0;
    _cxzW=0;_cyzW=0;
    _percent_of_hits_in_overlap_planes=0;

    bool Intrusion = false;
    int  entry_id = -1;
    std::vector<double> adcList;
    double ADCx = 0.;
    double ADCy = 0.;

    for (unsigned i = 0; i != slice.NCell(); ++i) {
      const art::Ptr<rb::CellHit> chit = slice.Cell(i);

      double tns = chit->TNS();
      int c = chit->Cell();
      int p = chit->Plane();
      adcList.emplace_back(chit->ADC());

      if(IsSurfaceHit(p,c)==true) _number_of_surface_hits++;
      _meanADC += chit->ADC();
      _cz += p;
      _czW += p * (chit->ADC());

      if (chit->View()==geo::kX) {       //X view
        Xhitlist.emplace_back(chit);
        ADCx += chit->ADC();
        ++_nxHits;
        _cx += c;
        _cxW += c*(chit->ADC());
        _cxz+=p;
        _cxzW+=p*(chit->ADC());

        if (PX_max < p) PX_max = p;
        if (PX_min > p) PX_min = p;
        if (CX_min > c) CX_min = c;
        if (CX_max < c) CX_max = c;
        if (TX_max < tns) TX_max = tns;
        if (TX_min > tns) TX_min = tns;

        if (!Intrusion) {
          if (c < _xdelt){
            entry_id  = 1;
            Intrusion = true;
          }
          else if (c > _xMax - _xdelt){
            entry_id  = 2;
            Intrusion = true;
          }
        }
        else if (!_penetrated) {
          if      (c < _xdelt && entry_id!=1)       _penetrated = true;
          else if (c > _xMax-_xdelt && entry_id!=2) _penetrated = true;
        }
      }

      else {       //Y view
        Yhitlist.emplace_back(chit);
        ADCy += chit->ADC();
        ++_nyHits;
        _cy += c;
        _cyW += c*(chit->ADC());
        _cyz+=p;
        _cyzW+=p*(chit->ADC());

        if (PY_max < p) PY_max = p;
        if (PY_min > p) PY_min = p;
        if (CY_min > c) CY_min = c;
        if (CY_max < c) CY_max = c;
        if (TY_max < tns) TY_max = tns;
        if (TY_min > tns) TY_min = tns;

        if (!Intrusion) {
          if (chit->Cell() < _ydelt) {
            entry_id  = 3;
            Intrusion = true;
          }
          else if (chit->Cell() > _yMax - _ydelt) {
            entry_id  = 4;
            Intrusion = true;
          }
        }

        else if (!_penetrated) {
          if      (chit->Cell() < _ydelt && entry_id!=3)       _penetrated = true;
          else if (chit->Cell() > _yMax-_ydelt && entry_id!=4) _penetrated = true;
        }
      }

      if (!Intrusion) {
        if (chit->Plane() < _zdelt) {
          entry_id  = 5;
          Intrusion = true;
        }
        else if (chit->Plane() > _zMax - _zdelt) {
          entry_id  = 6;
          Intrusion = true;
        }
      }
      else if (!_penetrated) {       //intrusion happened but not through yet
        if (chit->Plane() < _zdelt && entry_id!=5) _penetrated=true;
        else if (chit->Plane() > _zMax-_zdelt && entry_id!=6) _penetrated=true;
      }
    }
    //print out slice information:
    std::cout<<"X view Plane range: "<<PX_min<<" "<<PX_max<<std::endl;
    std::cout<<"X view Cell range: "<<CX_min<<" "<<CX_max<<std::endl;
    std::cout<<"X view time range: "<<TX_min<<" "<<TX_max<<std::endl;

    std::cout<<"Y view Plane range: "<<PY_min<<" "<<PY_max<<std::endl;
    std::cout<<"Y view Cell range: "<<CY_min<<" "<<CY_max<<std::endl;
    std::cout<<"Y view time range: "<<TY_min<<" "<<TY_max<<std::endl;


    //Record the info
    double nhits = _nxHits+_nyHits+0.;
    _cx  /= (_nxHits+0.);
    _cxW /= ADCx; 
    _cy  /= (_nyHits+0.);
    _cyW /= ADCy; 
    _cz  /= nhits;
    _czW /= _meanADC; 


    //    _meanADC   = _meanADC/nhits;
    _rmsADC    = RMS(adcList, _meanADC);

    _dpx          = PX_max-PX_min;
    _dpy          = PY_max-PY_min;
    _dcx          = CX_max-CX_min;
    _dcy          = CY_max-CY_min;
    _dp           =std::max(_dpx,_dpy);
    _tmin        = std::min(TX_min,TY_min);
    _tmax        = std::max(TX_max,TY_max);

    //preselection: if failed these, don't record.
    std::cout<<"nxHits: "<<_nxHits<<" nyHits: "<<_nyHits<<" penetration: "<<_penetrated<<std::endl;
    if(_nxHits<3||_nyHits<3||(!_penetrated)||(std::min(PX_max,PY_max) < std::max(PX_min,PY_min)+1)||(std::min(TX_max,TY_max)+62 < std::max(TX_min,TY_min) ))   continue;
    

    //added by enhao
    _cxz /= (_nxHits+0.0);
    _cxzW/=ADCx;
    _cyz /= (_nyHits+0.0);
    _cyzW/=ADCy;
    std::cout<<"number of surface hits: "<<_number_of_surface_hits<<std::endl;
    if(_number_of_surface_hits<22) continue;
    double _diagonal_length=std::sqrt(2.82072025*_dp*_dp+_dcx*_dcx+_dcy*_dcy);
    std::cout<<"slice diagonal length: "<<_diagonal_length<<std::endl;
    if(_diagonal_length<79) continue;
    _number_of_hits_per_length=1.0*(_nxHits+_nyHits)/std::sqrt(2.82072025*_dp*_dp+_dcx*_dcx+_dcy*_dcy);
    std::cout<<"number of hits divided by slice diagonal length"<<_number_of_hits_per_length<<std::endl;
    if(_number_of_hits_per_length<0.43) continue;
    double tracklength_of_xview=std::sqrt(2.82072025*_dpx*_dpx+_dcx*_dcx);
    double tracklength_of_yview=std::sqrt(2.82072025*_dpy*_dpy+_dcy*_dcy);
    NumberOfCellsPerLength(slice,tracklength_of_xview,tracklength_of_yview,_number_of_cells_per_length_xview,_number_of_cells_per_length_yview);
    std::cout<<"number of cells per length xview: "<<_number_of_cells_per_length_xview<<" yview: "<<_number_of_cells_per_length_yview<<std::endl;
    if(_number_of_cells_per_length_xview<0.19||_number_of_cells_per_length_yview<0.19)continue;
    StdevCellsPerPlane(slice, _stdev_cells_per_plane_xview, _stdev_cells_per_plane_yview);
    std::cout<<"stdev cells per plane xview: "<<_stdev_cells_per_plane_xview<<" yview: "<<_stdev_cells_per_plane_yview<<std::endl;
    if(_stdev_cells_per_plane_xview>1.26||_stdev_cells_per_plane_yview>1.21)continue;
    _weighted_off_center_xx=std::abs((_cxW-0.5*(CX_max+CX_min))/_dcx); 
    _weighted_off_center_xz=std::abs((_cxzW-0.5*(PX_max+PX_min))/_dpx); 
    _weighted_off_center_yy=std::abs((_cyW-0.5*(CY_max+CY_min))/_dcy); 
    _weighted_off_center_yz=std::abs((_cyzW-0.5*(PY_max+PY_min))/_dpy); 
    _percent_of_hits_in_overlap_planes=NumberOfHitsInOverlapPlanesCut(slice,PX_min,PX_max,PY_min,PY_max);


    std::cout<<"percent of hits in overlap planes: "<<_percent_of_hits_in_overlap_planes<<std::endl;
    std::cout<<"weighted off center xz: "<<_weighted_off_center_xz<<" yz: "<<_weighted_off_center_yz<<std::endl;
    if(_percent_of_hits_in_overlap_planes>0.45&&_weighted_off_center_xz<0.15&&_weighted_off_center_yz<0.15){
      std::cout<<"trigger!!!!!!!!!!!!!"<<std::endl;
      std::cout<<"trigger slice diagonal length:  "<<std::sqrt(2.82072025*_dp*_dp+_dcx*_dcx+_dcy*_dcy)<<std::endl;
      std::cout<<"evt id: "<<evtID<<std::endl;
    clustercol->push_back(slice);
    }
  }
    evt.put(std::move(clustercol));

} //end analyzing 

double zcl::FmmTriggerEvd::RMS(std::vector<double> & v, double & mean) {
  double rms = 0;
  mean = 0;

  for (unsigned i = 0; i!= v.size(); ++i) {
    mean += v[i];
  }
  mean /= double(v.size());
  for (unsigned i = 0; i!= v.size(); ++i) {
    rms += (v[i]-mean)*(v[i]-mean);
  }
  rms /= double(v.size());
  rms = sqrt(rms);
  return rms;
}

bool zcl::FmmTriggerEvd::IsSurfaceHit(int plane,int cell) const{
  if (plane<_zdelt) return true;
  if (plane>_zMax-_zdelt)return true;
  if (cell<_xdelt) return true;
  if (cell>_xMax-_xdelt) return true;
  return false;
}


void zcl::FmmTriggerEvd::NumberOfCellsPerLength(rb::Cluster const &slice,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 (unsigned i = 0; i != slice.NCell(); ++i) {
    const art::Ptr<rb::CellHit> chit = slice.Cell(i);
    int c = chit->Cell();
    int p = chit->Plane();
    if (chit->View()==geo::kX) {
      x_hits_plane[p-PX_min].push_back(c); 
    }
    else {
      y_hits_plane[p-PY_min].push_back(c); 
    }
  }
  //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 (size_t 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 (size_t 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 (size_t 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 (size_t 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 zcl::FmmTriggerEvd::StdevCellsPerPlane(rb::Cluster const &slice, double &stdev_cells_per_plane_xview, double &stdev_cells_per_plane_yview ) const{
  //fluctuation of number of cells of each plane check by Enhao
  //I need to use cells not hits, because there might be many hits in same cell, for example sag effect, and the cosmic shower has many different cells hits
  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 (unsigned i = 0; i != slice.NCell(); ++i) {
    const art::Ptr<rb::CellHit> chit = slice.Cell(i);
    int c = chit->Cell();
    int p = chit->Plane();
    if (chit->View()==geo::kX) {
      x_hits_plane[p-PX_min].push_back(c); 
    }
    else {
      y_hits_plane[p-PY_min].push_back(c); 
    }
  }
  double sum =0;
  double mean = 0;
  double sq_sum = 0;
  double 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 (size_t 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 (size_t 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;
}
double zcl::FmmTriggerEvd::NumberOfHitsInOverlapPlanesCut(rb::Cluster const &slice,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 (unsigned i = 0; i != slice.NCell(); ++i) {
    const art::Ptr<rb::CellHit> chit = slice.Cell(i);
    if(chit->Plane()<=max && chit->Plane()>=min)    number_hits_in_overlap_planes++;
  }
  return number_hits_in_overlap_planes/slice.NCell();
}
DEFINE_ART_MODULE(zcl::FmmTriggerEvd)