FastMonopoleTriggers.cxx

Go to the documentation of this file.

#include "FastMonopoleTriggers.h"

bool zcl::FastMonopoleTriggers::is_trigger_by_epoch1_fmmtrigger(art::Handle< std::vector<rb::Cluster> > slices, int & _penetrated, double & _nSatHits, double & _meanADC){
  bool triggered = false;
  for (unsigned sliceIdx = 0; sliceIdx != slices->size(); ++sliceIdx) {
    const rb::Cluster& slice = (*slices)[sliceIdx];
    double _nxhits_slice = 0;
    double _nyhits_slice = 0;
    double _deltaP;
    double _deltaTNS;
    double PX_max = 0;  
    double PY_max = 0;
    double PX_min = 999;
    double PY_min = 999;
    double TX_max = -9e7;       
    double TY_max = -9e7;
    double TX_min = 9e9;
    double TY_min = 9e9;

    _penetrated = 0;
    double _trueHits   = 0;
    _nSatHits   = 0;
    _meanADC    = 0;
    //double noise_hits  = 0;

    bool Intrusion = false;
    int entry_id = -1;
    for (unsigned i = 0; i != slice.NCell(); ++i) {
      const art::Ptr<rb::CellHit> chit = slice.Cell(i);
      if(chit->IsMC()==true) ++_trueHits;
    }

    if(_trueHits==0) continue;
    for (unsigned i = 0; i != slice.NCell(); ++i) {
      const art::Ptr<rb::CellHit> chit = slice.Cell(i);
      if (chit->ADC() < 501) continue;
      if (chit->ADC() > 3100) ++_nSatHits;
      _meanADC += chit->ADC();
      float tns = chit->TNS();

      if (chit->View()==geo::kX) {       //X view
        ++_nxhits_slice;
        if (PX_max < chit->Plane() ) PX_max = chit->Plane();
        if (PX_min > chit->Plane() ) PX_min = chit->Plane();
        if (TX_max < tns   ) TX_max = tns;
        if (TX_min > tns   ) TX_min = tns;

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

      else {       //Y view
        ++_nyhits_slice;
        if (PY_max < chit->Plane() ) PY_max = chit->Plane();
        if (PY_min > chit->Plane() ) PY_min = chit->Plane();
        if (TY_max < chit->TNS()   ) TY_max = chit->TNS();
        if (TY_min > chit->TNS()   ) TY_min = chit->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 = 1;
          else if (chit->Cell() > _yMax-_ydelt && entry_id!=4) _penetrated = 1;
        }
      }

      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=1;
        else if (chit->Plane() > _zMax-_zdelt && entry_id!=6) _penetrated=1;
      }
    }
    //Record the info
    _meanADC     = _meanADC/(_nxhits_slice+_nyhits_slice+0.);
    _deltaP      = std::min(PX_max,PY_max) - std::max(PX_min,PY_min);
    _deltaTNS    = std::min(TX_max,TY_max) - std::max(TX_min,TY_min);    

    if(_trueHits>2&&_penetrated==1&&_nSatHits>4&&_meanADC>1200&&_deltaP>2&&_deltaTNS>0&&_nxhits_slice>2&&_nyhits_slice>2){
      triggered = true;
      return triggered;
    }
  }
  _penetrated = 0;
  _nSatHits   = 0;
  _meanADC    = 0;
  return triggered;
} 


bool zcl::FastMonopoleTriggers::passed_epoch2_fmmtrigger(art::Handle< std::vector<rb::Cluster> > slices, int & _penetrated, double & _nSatHits, double & _meanADC){
  bool triggered = false;
  for (unsigned sliceIdx = 0; sliceIdx != slices->size(); ++sliceIdx) {
    const rb::Cluster& slice = (*slices)[sliceIdx];
    int nx = 0;
    int ny = 0;
    double sumADC=0;
    //    double stdevADC=0;
    double sq_sum_ADC=0;
    double _deltaP;
    double _deltaTNS;
    double PX_max = 0;  
    double PY_max = 0;
    double PX_min = 9999;
    double PY_min = 9999;
    double CX_max = 0;  
    double CY_max = 0;
    double CX_min = 9999;
    double CY_min = 9999;
    double TX_max = -9e7;       
    double TY_max = -9e7;
    double TX_min = 9e9;
    double TY_min = 9e9;
    double _trueHits   = 0;

    _penetrated = 0;
    _nSatHits   = 0;
    _meanADC    = 0;
    //double noise_hits  = 0;

    bool Intrusion = false;
    int entry_id = -1;
    for (unsigned i = 0; i != slice.NCell(); ++i) {
      const art::Ptr<rb::CellHit> chit = slice.Cell(i);
      if(chit->IsMC()==true) ++_trueHits;
    }

    if(_trueHits==0) continue;
    for (unsigned i = 0; i != slice.NCell(); ++i) {
      const art::Ptr<rb::CellHit> chit = slice.Cell(i);
      if (chit->ADC() < 501) continue;
      if (chit->ADC() > 3100) ++_nSatHits;
      sumADC += chit->ADC();
      sq_sum_ADC += chit->ADC()*chit->ADC();
      _meanADC += chit->ADC();
      float tns = chit->TNS();

      if (chit->View()==geo::kX) {       //X view
        ++nx;
        if (PX_max < chit->Plane() ) PX_max = chit->Plane();
        if (PX_min > chit->Plane() ) PX_min = chit->Plane();
        if (CX_max < chit->Cell() )  CX_max  = chit->Cell();
        if (CX_min > chit->Cell() )  CX_min  = chit->Cell();
        if (TX_max < tns   ) TX_max = tns;
        if (TX_min > tns   ) TX_min = tns;

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

      else {       //Y view
        ++ny;
        if (PY_max < chit->Plane() ) PY_max = chit->Plane();
        if (PY_min > chit->Plane() ) PY_min = chit->Plane();
        if (CY_max < chit->Cell()  ) CY_max = chit->Cell();
        if (CY_min > chit->Cell()  ) CY_min = chit->Cell();
        if (TY_max < chit->TNS()   ) TY_max = chit->TNS();
        if (TY_min > chit->TNS()   ) TY_min = chit->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 = 1;
          else if (chit->Cell() > _yMax-_ydelt && entry_id!=4) _penetrated = 1;
        }
      }

      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=1;
        else if (chit->Plane() > _zMax-_zdelt && entry_id!=6) _penetrated=1;
      }
    }
    //Record the info
    _meanADC     = _meanADC/(nx+ny+0.);
    _deltaP      = std::min(PX_max,PY_max) - std::max(PX_min,PY_min);
    _deltaTNS    = std::min(TX_max,TY_max) - std::max(TX_min,TY_min);    

    if(_trueHits>2&&_penetrated==1&&_deltaP>2&&_deltaTNS>0&&nx>2&&ny>2){
      float tracklength =Distance(CX_max, CY_max, PX_max, CX_min, CY_min, PX_min); 
      if( tracklength<80) continue;
      if( NumberOfSurfaceHits(slice) < 22) continue; 
      if( slice.NCell()*1.0/tracklength < 0.44) continue;
      
      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);
      double number_of_cells_per_length_xview=0;
      double number_of_cells_per_length_yview=0;
      double weighted_off_center_xx=0, weighted_off_center_xz=0, weighted_off_center_yy=0, weighted_off_center_yz=0;
      double stdev_cells_per_plane_xview=0, stdev_cells_per_plane_yview=0;

      NumberOfCellsPerLength(slice, tracklength_of_xview,tracklength_of_yview,number_of_cells_per_length_xview,number_of_cells_per_length_yview, PX_max, PY_max, PX_min, PY_min);
      if( number_of_cells_per_length_xview<0.2 || number_of_cells_per_length_yview<0.2) continue;

      //stdevADC = std::sqrt( sq_sum_ADC*1.0/(nx+ny) - sumADC*1.0/(nx+ny)*sumADC*1.0/(nx+ny));
      StdevCellsPerPlane(slice, stdev_cells_per_plane_xview, stdev_cells_per_plane_yview, PX_max, PX_min, PY_max, PY_min );
      if (stdev_cells_per_plane_xview > 1.25 || stdev_cells_per_plane_yview > 1.2) continue;
      WeightedCenterCut(slice, 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 (weighted_off_center_xz > 0.13 || weighted_off_center_yz > 0.13) continue; 
      //sumADC*1.0/(nx+ny);
      //stdevADC;
      //nx;
      //ny;
      if( NumberOfHitsInOverlapPlanesCut(slice,PX_min,PX_max,PY_min,PY_max) <= 0.45) continue;
      triggered = true;
      return triggered;
    }
  }
  _penetrated = 0;
  _nSatHits   = 0;
  _meanADC    = 0;
  return triggered;
} 

void zcl::FastMonopoleTriggers::NumberOfCellsPerLength(rb::Cluster const &hits,double tracklength_of_xview, double tracklength_of_yview, double &number_of_cells_per_length_xview,double &number_of_cells_per_length_yview, double const &PX_max, double const  &PY_max, double const &PX_min, double const  &PY_min) const{
  //Maybe use hits instead of cells?
  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_ptr: hits.AllCells()) {
    const rb::CellHit hit = *hit_ptr;
    if (hit.View() == geo::kX) {
      //std::cout<<"X hits"<<std::endl;
      x_hits_plane[hit.Plane()-PX_min].push_back(hit.Cell()); 
    }
    if (hit.View() == geo::kY) {
      //std::cout<<"Y hits"<<std::endl;
      y_hits_plane[hit.Plane()-PY_min].push_back(hit.Cell()); 
    }
  }
  //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;

}

int zcl::FastMonopoleTriggers::NumberOfSurfaceHits(rb::Cluster const &hits) const{
  int number_of_surface_hits=0;
  for(unsigned i=0;i!=hits.NCell();++i){
    if(SurfAssign(*(hits.Cell(i)))==true) number_of_surface_hits++;
  }
  return number_of_surface_hits;
}

bool zcl::FastMonopoleTriggers::SurfAssign(rb::CellHit const &hit) const{
  //Check if this is close to surface
  if      ( hit.Cell()  > _xMax - _xdelt   ) return true;
  else if ( hit.Cell()  < _xdelt     ) return true;
  else if ( hit.Plane() < _zdelt     ) return true;
  else if ( hit.Plane() > _zMax - _zdelt   ) return true;

  return false;
}
void zcl::FastMonopoleTriggers::StdevCellsPerPlane(rb::Cluster const &hits, double &stdev_cells_per_plane_xview, double &stdev_cells_per_plane_yview, double PX_max, double PX_min, double PY_max, double PY_min ) 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, and the cosmic shower has many different cells hits
  //Maybe use hits instead of cells?
  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_ptr: hits.AllCells()) {
    const rb::CellHit hit = *hit_ptr;
    if (hit.View() == geo::kX) {
      //std::cout<<"X hits"<<std::endl;
      x_hits_plane[hit.Plane()-PX_min].push_back(hit.Cell()); 
    }
    if (hit.View() == geo::kY) {
      //std::cout<<"Y hits"<<std::endl;
      y_hits_plane[hit.Plane()-PY_min].push_back(hit.Cell()); 
    }
  }
  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 (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);
  //        std::cout<<"x stdev is "<<stdev<<" "<<stdev/mean<<std::endl;
  //if (stdev/mean > _stdevcellsperplane) return false;
  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;
  //        std::cout<<"y stdev is "<<stdev<<" "<<stdev/mean<<std::endl;
  //        std::cout<<"sum of x and y stdev is "<<sum_stdev<<std::endl;
  //if (stdev/mean > _stdevcellsperplane) return false;
}

float zcl::FastMonopoleTriggers::Distance( double const &xcell1, double const  &ycell1, double const &plane1,double const &xcell2, double const &ycell2 , double const &plane2) const {
  return std::sqrt((plane1-plane2)*(plane1-plane2)*2.82072+(xcell1-xcell2)*(xcell1-xcell2)+(ycell1-ycell2)*(ycell1-ycell2));
  //This part I should use Martin's smt::Constant::cell width, plane width
}

double zcl::FastMonopoleTriggers::WeightedCenterCut(rb::Cluster 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_ptr: hits.AllCells()) {
    const rb::CellHit hit = *hit_ptr;
    if (hit.View() == geo::kX) {
      SumOfXHitsADC+=hit.ADC();
      Weighted_PX_Center+=hit.ADC()*(hit.Plane()-PX_min);
      Weighted_CX_Center+=hit.ADC()*(hit.Cell()-CX_min);
    }
    if (hit.View() == geo::kY) {
      SumOfYHitsADC+=hit.ADC();
      Weighted_PY_Center+=hit.ADC()*(hit.Plane()-PY_min);
      Weighted_CY_Center+=hit.ADC()*(hit.Cell()-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);
  //  std::cout<<"weighted px, py, cx, cy center is "<<Weighted_PX_Center<<" "<<Weighted_PY_Center<<" "<<Weighted_CX_Center<<" "<<Weighted_CY_Center<<std::endl;
  // if(std::abs(Weighted_PX_Center-0.5)>_MaxOffCenter||std::abs(Weighted_PY_Center-0.5)>_MaxOffCenter||std::abs(Weighted_CX_Center-0.5)>_MaxOffCenter||std::abs(Weighted_CY_Center-0.5)>_MaxOffCenter) return false;
  //else  return true;
  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 zcl::FastMonopoleTriggers::NumberOfHitsInOverlapPlanesCut(rb::Cluster 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_ptr: hits.AllCells()) {
    const rb::CellHit hit = *hit_ptr;
    if(hit.Plane()<=max && hit.Plane()>=min)    number_hits_in_overlap_planes++;
  }

  return number_hits_in_overlap_planes/hits.NCell();

}