FastMMTrigger_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // Class: FastMMTrigger
3 // Module Type: filter
4 // File: FastMMTrigger_module.cc
5 //
6 // Created at Nov 22 2015 by Enhao Song
7 ////////////////////////////////////////////////////////////////////////
18 #include <fhiclcpp/ParameterSet.h>
21 
30 
31 #include <string>
32 #include <vector>
33 #include <algorithm>
34 #include <numeric>
35 #include <cmath>
36 #include <map>
37 
38 namespace novaddt {
39  class FastMMTrigger;
40 }
42  public:
43  explicit FastMMTrigger(fhicl::ParameterSet const & p);
44  virtual ~FastMMTrigger();
45  virtual bool filter(art::Event & e) override;
46  void StdevCellsPerPlane(HitList const &hits, double &stdev_cells_per_plane_xview, double &stdev_cells_per_plane_yview ) const;
47  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;
48  float Distance(unsigned const &xcell1, unsigned const &ycell1, unsigned const &plane1,unsigned const &xcell2, unsigned const &ycell2 , unsigned const &plane2) const;
49  bool SurfAssign(DAQHit const & hit) const;
50  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;
51  double NumberOfHitsInOverlapPlanesCut(HitList const &hits,float PX_min,float PX_max,float PY_min,float PY_max ) const;
52  int NumberOfSurfaceHits(HitList const &hits) const;
53 
54  private:
55  unsigned _prescale;
58  int _xMin;
59  int _xMax;
60  int _yMin;
61  int _yMax;
62  int _zMin;
63  int _zMax;
64  int _deltx;
65  int _delty;
66  int _deltz;
67  int _sigmaP;
69 
72  unsigned _trigger_counts;
73 };
75  :
76  _prescale(p.get<unsigned>("prescale")),
77  _slicelabel(p.get<std::string>("slice_label")),
78  _sliceinstance(p.get<std::string>("slice_instance")),
79  _xMin(p.get<int>("x_min")),
80  _xMax(p.get<int>("x_max")),
81  _yMin(p.get<int>("y_min")),
82  _yMax(p.get<int>("y_max")),
83  _zMin(p.get<int>("z_min")),
84  _zMax(p.get<int>("z_max")),
85  _deltx(p.get<int>("x_delt")),
86  _delty(p.get<int>("y_delt")),
87  _deltz(p.get<int>("z_delt")),
88  _sigmaP(p.get<int>("sigma_P")),
89  _Min_number_of_hits_per_length(p.get<double>("Min_number_of_hits_per_length")),
90  _Min_number_of_cells_per_length_xview(p.get<double>("Min_number_of_cells_per_length_xview")),
91  _Min_number_of_cells_per_length_yview(p.get<double>("Min_number_of_cells_per_length_yview")),
92  _Max_stdev_cells_per_plane_xview(p.get<double>("Max_stdev_cells_per_plane_xview")),
93  _Max_stdev_cells_per_plane_yview(p.get<double>("Max_stdev_cells_per_plane_yview")),
94  _Max_weighted_off_center_yz(p.get<double>("Max_weighted_off_center_yz")),
95  _Max_weighted_off_center_xz(p.get<double>("Max_weighted_off_center_xz")),
96  _Min_percent_of_hits_in_overlap_planes(p.get<double>("Min_percent_of_hits_in_overlap_planes")),
97  _Min_number_of_surface_hits(p.get<double>("Min_number_of_surface_hits")),
98  _Min_slice_diagonal_length(p.get<double>("Min_slice_diagonal_length")),
100 {
101  // Call appropriate Produces<>() functions here.
102  produces<std::vector<TriggerDecision>>();
103 }
104 
106  // Clean up dynamic memory and other resources here.
107 }
108 
110  bool result = false;
111  std::unique_ptr<std::vector<TriggerDecision>>
112  trigger_decisions(new std::vector<TriggerDecision>);
115 
116  for (unsigned i = 0; i!= slices->size(); ++i) {
117  HitList const &hits = (*slices)[i];
118  PX_max = 0;//P is for plane
119  CX_max = 0;
120  CY_max = 0;
121  PY_max = 0;
122  PX_min = 99999;
123  CX_min = 99999;
124  CY_min = 99999;
125  PY_min = 99999;
126  TX_max = hits.front().TDC().val;
127  TY_max = TX_max;
128  TX_min = hits.back().TDC().val;
129  TY_min = TX_min;
130 
131  //int nb_hits = hits.size();
132  //std::cout<<"slice number = "<<i+1<<" number of hits: "<<nb_hits<<std::endl;
133  int entry_id = 0;
134  bool PassThrough = false;
135  bool Intrusion = false;
136 
137  int nx = 0;
138  int ny = 0;
139 
140  for (auto const hit: hits) {
141  // std::cout<<"TDC: "<<hit.TDC().val << " ADC: "<<hit.ADC().val<<" Plane: "<<hit.Plane().val<<" Cell: "<<hit.Cell().val<<std::endl;
142  //X view hits
143  if (hit.View().val == daqchannelmap::X_VIEW) {
144  //std::cout<<"X hits"<<std::endl;
145  ++nx;
146  if (PX_max < hit.Plane().val) PX_max = hit.Plane().val;
147  if (PX_min > hit.Plane().val) PX_min = hit.Plane().val;
148  if (CX_max < hit.Cell().val) CX_max = hit.Cell().val;
149  if (CX_min > hit.Cell().val) CX_min = hit.Cell().val;
150  if (TX_max < hit.TDC().val) TX_max = hit.TDC().val;
151  if (TX_min > hit.TDC().val) TX_min = hit.TDC().val;
152 
153  if (!Intrusion) {
154  if (hit.Cell().val < _xMin+_deltx) {
155  entry_id = 1;
156  Intrusion = true;
157  }
158  else if (hit.Cell().val> _xMax-_deltx) {
159  entry_id = 2;
160  Intrusion = true;
161  }
162  }
163  else if (!PassThrough) {
164  if (hit.Cell().val < _xMin+_deltx && entry_id != 1) PassThrough = true;
165  else if (hit.Cell().val> _xMax-_deltx && entry_id != 2) PassThrough = true;
166  }
167  }
168  //Y view hits
169  else {
170  //std::cout<<"Y hits"<<std::endl;
171  ++ny;
172  if (PY_max < hit.Plane().val) PY_max = hit.Plane().val;
173  if (PY_min > hit.Plane().val) PY_min = hit.Plane().val;
174  if (CY_max < hit.Cell().val) CY_max = hit.Cell().val;
175  if (CY_min > hit.Cell().val) CY_min = hit.Cell().val;
176  if (TY_max < hit.TDC().val ) TY_max = hit.TDC().val;
177  if (TY_min > hit.TDC().val ) TY_min = hit.TDC().val;
178 
179  if (!Intrusion) {
180  if (hit.Cell().val < _yMin+_delty) {
181  entry_id = 3;
182  Intrusion = true;
183  }
184  else if (hit.Cell().val> _yMax-_delty) {
185  entry_id = 4;
186  Intrusion = true;
187  }
188  }
189  else if (!PassThrough) {
190  if (hit.Cell().val < _yMin+_delty && entry_id != 3) PassThrough = true;
191  else if (hit.Cell().val> _yMax-_delty && entry_id != 4) PassThrough = true;
192  }
193  }
194  //No matter which view:
195  if (!Intrusion) {
196  if (hit.Plane().val < _zMin+_deltz) {
197  entry_id = 5;
198  Intrusion = true;
199  }
200  else if (hit.Plane().val > _zMax-_deltz) {
201  entry_id = 6;
202  Intrusion = true;
203  }
204  }
205  else if (!PassThrough) {
206  if (hit.Plane().val < _zMin+_deltz && entry_id != 5) PassThrough=true;
207  else if (hit.Plane().val > _zMax-_deltz && entry_id != 6) PassThrough=true;
208  }
209  }
210 
211  //Min-hits on both view:
212  // std::cout<<"xhits number : "<<nx<<" yhits number : "<<ny<<std::endl;
213  if (nx<10 || ny<10) continue;
214  // std::cout<<"minimum hits requirement met!"<<std::endl;
215  //Passing Through Test:
216  if (!PassThrough) continue;
217  // std::cout<<"through!"<<std::endl;
218  //std::cout<<"X view Plane range: "<<PX_min<<" "<<PX_max<<std::endl;
219  //std::cout<<"X view Cell range: "<<CX_min<<" "<<CX_max<<std::endl;
220  //std::cout<<"X view time range: "<<TX_min<<" "<<TX_max<<std::endl;
221 
222  //std::cout<<"Y view Plane range: "<<PY_min<<" "<<PY_max<<std::endl;
223  //std::cout<<"Y view Cell range: "<<CY_min<<" "<<CY_max<<std::endl;
224  //std::cout<<"Y view time range: "<<TY_min<<" "<<TY_max<<std::endl;
225 
226  //Matching Test:
227  if (std::min(PX_max,PY_max) < std::max(PX_min,PY_min)+_sigmaP) continue; //Require minimum overlapping planes
228  if (std::min(TX_max,TY_max) < std::max(TX_min,TY_min) ) continue;
229 
230  unsigned dT = std::max(TX_max,TY_max)-std::min(TX_min,TY_min);
231 
232  float tracklength =Distance(CX_max, CY_max, PX_max, CX_min, CY_min, PX_min);
233  //std::cout<<"number of surface hits: "<<NumberOfSurfaceHits(hits)<<std::endl;
235  //std::cout<<"slice diagonal length: "<<tracklength<<std::endl;
236  if(tracklength<_Min_slice_diagonal_length) continue;
237  //std::cout<<"number of hits divided by slice diagonal length"<<hits.size()/tracklength<<std::endl;
238  if(_Min_number_of_hits_per_length>hits.size()/tracklength) continue;
239 
240  double number_of_cells_per_length_xview=0;
241  double number_of_cells_per_length_yview=0;
242  double tracklength_of_xview=Distance(CX_max, 0, PX_max, CX_min, 0, PX_min);
243  double tracklength_of_yview=Distance(0, CY_max, PY_max, 0, CY_min, PY_min);
244  NumberOfCellsPerLength(hits,tracklength_of_xview,tracklength_of_yview,number_of_cells_per_length_xview,number_of_cells_per_length_yview);
245  //std::cout<<"number of cells per length xview: "<<number_of_cells_per_length_xview<<" yview: "<<number_of_cells_per_length_yview<<std::endl;
246  if(_Min_number_of_cells_per_length_xview>number_of_cells_per_length_xview) continue;
247  if(_Min_number_of_cells_per_length_yview>number_of_cells_per_length_yview) continue;
248 
249  double stdev_cells_per_plane_xview=0, stdev_cells_per_plane_yview=0;
250  StdevCellsPerPlane(hits, stdev_cells_per_plane_xview, stdev_cells_per_plane_yview);
251  //std::cout<<"stdev cells per plane xview: "<<stdev_cells_per_plane_xview<<" yview: "<<stdev_cells_per_plane_yview<<std::endl;
252  if(_Max_stdev_cells_per_plane_xview<stdev_cells_per_plane_xview) continue;
253  if(_Max_stdev_cells_per_plane_yview<stdev_cells_per_plane_yview) continue;
254 
255  double weighted_off_center_xx=0, weighted_off_center_xz=0, weighted_off_center_yy=0, weighted_off_center_yz=0;
256  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 );
257  //std::cout<<"weighted off center xz: "<<weighted_off_center_xz<<" yz: "<<weighted_off_center_yz<<std::endl;
258  if(_Max_weighted_off_center_xz<=weighted_off_center_xz) continue;
259  if(_Max_weighted_off_center_yz<=weighted_off_center_yz) continue;
260 
261  //std::cout<<"percent of hits in overlap planes: "<<NumberOfHitsInOverlapPlanesCut(hits,PX_min,PX_max,PY_min,PY_max)<<std::endl;
263 
264  ++_trigger_counts;
265 
267  trigger_decisions->
268  emplace_back(std::min(TX_min,TY_min),
269  dT,
271  result = true;
272  }
273  }
274  e.put(std::move(trigger_decisions));
275 
276  return result;
277 }
278 
279 
280 void novaddt::FastMMTrigger::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{
281  std::vector<std::vector<int>> x_hits_plane(PX_max-PX_min+1,std::vector<int>(0));
282  std::vector<std::vector<int>> y_hits_plane(PY_max-PY_min+1,std::vector<int>(0));
283  for (auto const hit: hits) {
284  if (hit.View().val == daqchannelmap::X_VIEW) {
285  //std::cout<<"X hits"<<std::endl;
286  x_hits_plane[hit.Plane().val-PX_min].push_back(hit.Cell().val);
287  }
288  if (hit.View().val == daqchannelmap::Y_VIEW) {
289  //std::cout<<"Y hits"<<std::endl;
290  y_hits_plane[hit.Plane().val-PY_min].push_back(hit.Cell().val);
291  }
292  }
293  //following is to get unique number of cells hit per plane
294  std::vector<int> x_cells_plane(PX_max-PX_min+1,0);
295  std::vector<int> y_cells_plane(PY_max-PY_min+1,0);
296  for (unsigned i=0;i!=x_hits_plane.size();i++){
297  std::sort(x_hits_plane[i].begin(),x_hits_plane[i].end());
298  std::vector<int>::iterator it;
299  it = std::unique (x_hits_plane[i].begin(),x_hits_plane[i].end());
300  x_hits_plane[i].resize(std::distance(x_hits_plane[i].begin(),it));
301  x_cells_plane[i]=x_hits_plane[i].size();
302  }
303  for (unsigned i=0;i!=y_hits_plane.size();i++){
304  std::sort(y_hits_plane[i].begin(),y_hits_plane[i].end());
305  std::vector<int>::iterator it;
306  it = std::unique (y_hits_plane[i].begin(),y_hits_plane[i].end());
307  y_hits_plane[i].resize(std::distance(y_hits_plane[i].begin(),it));
308  y_cells_plane[i]=y_hits_plane[i].size();
309  }
310  //following is to get number of cells per track length in x view and y view
311  int sum_of_cells=0;
312  for (unsigned i=0;i!=x_cells_plane.size();i++){
313  sum_of_cells+=x_cells_plane[i];
314  }
315  number_of_cells_per_length_xview=sum_of_cells/tracklength_of_xview;
316  sum_of_cells=0;
317  for (unsigned i=0;i!=y_cells_plane.size();i++){
318  sum_of_cells+=y_cells_plane[i];
319  }
320  number_of_cells_per_length_yview=sum_of_cells/tracklength_of_yview;
321 
322  }
323 void novaddt::FastMMTrigger::StdevCellsPerPlane(HitList const &hits, double &stdev_cells_per_plane_xview, double &stdev_cells_per_plane_yview ) const{
324  //fluctuation of number of cells of each plane check by Enhao
325  std::vector<std::vector<int>> x_hits_plane(PX_max-PX_min+1,std::vector<int>(0));
326  std::vector<std::vector<int>> y_hits_plane(PY_max-PY_min+1,std::vector<int>(0));
327  for (auto const hit: hits) {
328  if (hit.View().val == daqchannelmap::X_VIEW) {
329  //std::cout<<"X hits"<<std::endl;
330  x_hits_plane[hit.Plane().val-PX_min].push_back(hit.Cell().val);
331  }
332  if (hit.View().val == daqchannelmap::Y_VIEW) {
333  //std::cout<<"Y hits"<<std::endl;
334  y_hits_plane[hit.Plane().val-PY_min].push_back(hit.Cell().val);
335  }
336  }
337  double sum =0;
338  double mean = 0;
339  double sq_sum = 0;
340  double stdev = 0;
341  //double min_stdev=0;
342  //following is to get unique number of cells hit per plane
343  std::vector<int> x_cells_plane(PX_max-PX_min+1,0);
344  std::vector<int> y_cells_plane(PY_max-PY_min+1,0);
345  for (unsigned i=0;i!=x_hits_plane.size();i++){
346  std::sort(x_hits_plane[i].begin(),x_hits_plane[i].end());
347  std::vector<int>::iterator it;
348  it = std::unique (x_hits_plane[i].begin(),x_hits_plane[i].end());
349  x_hits_plane[i].resize(std::distance(x_hits_plane[i].begin(),it));
350  x_cells_plane[i]=x_hits_plane[i].size();
351  }
352  for (unsigned i=0;i!=y_hits_plane.size();i++){
353  std::sort(y_hits_plane[i].begin(),y_hits_plane[i].end());
354  std::vector<int>::iterator it;
355  it = std::unique (y_hits_plane[i].begin(),y_hits_plane[i].end());
356  y_hits_plane[i].resize(std::distance(y_hits_plane[i].begin(),it));
357  y_cells_plane[i]=y_hits_plane[i].size();
358  }
359  sum = std::accumulate(x_cells_plane.begin(),x_cells_plane.end(),0.0);
360  mean = sum / x_cells_plane.size();
361  sq_sum = std::inner_product(x_cells_plane.begin(),x_cells_plane.end(),x_cells_plane.begin(),0.0);
362  stdev = std::sqrt(sq_sum / x_cells_plane.size()-mean * mean);
363  stdev_cells_per_plane_xview=stdev/mean;
364  sum = std::accumulate(y_cells_plane.begin(),y_cells_plane.end(),0.0);
365  mean = sum / y_cells_plane.size();
366  sq_sum = std::inner_product(y_cells_plane.begin(),y_cells_plane.end(),y_cells_plane.begin(),0.0);
367  stdev = std::sqrt(sq_sum / y_cells_plane.size()-mean * mean);
368  stdev_cells_per_plane_yview=stdev/mean;
369 }
370 float novaddt::FastMMTrigger::Distance( unsigned const &xcell1, unsigned const &ycell1, unsigned const &plane1,unsigned const &xcell2, unsigned const &ycell2 , unsigned const &plane2) const {
371  return std::sqrt((plane1-plane2)*(plane1-plane2)*2.82072+(xcell1-xcell2)*(xcell1-xcell2)+(ycell1-ycell2)*(ycell1-ycell2));
372 }
373 
375  //Check if this is close to surface
376  if ( hit.Cell().val > _xMax - _deltx ) return true;
377  else if ( hit.Cell().val < _deltx ) return true;
378  else if ( hit.Plane().val < _deltz ) return true;
379  else if ( hit.Plane().val > _zMax - _deltz ) return true;
380 
381  return false;
382 }
384  int number_of_surface_hits=0;
385  for(unsigned i=0;i!=hits.size();++i){
386  if(SurfAssign(hits[i])==true) number_of_surface_hits++;
387  }
388  return number_of_surface_hits;
389 }
390 double novaddt::FastMMTrigger::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{
391  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;
392  for (auto const hit: hits) {
393  if (hit.View().val == daqchannelmap::X_VIEW) {
394  SumOfXHitsADC+=hit.ADC().val;
395  Weighted_PX_Center+=hit.ADC().val*(hit.Plane().val-PX_min);
396  Weighted_CX_Center+=hit.ADC().val*(hit.Cell().val-CX_min);
397  }
398  if (hit.View().val == daqchannelmap::Y_VIEW) {
399  SumOfYHitsADC+=hit.ADC().val;
400  Weighted_PY_Center+=hit.ADC().val*(hit.Plane().val-PY_min);
401  Weighted_CY_Center+=hit.ADC().val*(hit.Cell().val-CY_min);
402  }
403  }
404  if(PX_max!=PX_min) Weighted_PX_Center=Weighted_PX_Center/SumOfXHitsADC/(PX_max-PX_min);
405  if(PY_max!=PY_min) Weighted_PY_Center=Weighted_PY_Center/SumOfYHitsADC/(PY_max-PY_min);
406  if(CX_max!=CX_min) Weighted_CX_Center=Weighted_CX_Center/SumOfXHitsADC/(CX_max-CX_min);
407  if(CY_max!=CY_min) Weighted_CY_Center=Weighted_CY_Center/SumOfYHitsADC/(CY_max-CY_min);
408  weighted_off_center_xx=std::abs(Weighted_CX_Center-0.5);
409  weighted_off_center_xz=std::abs(Weighted_PX_Center-0.5);
410  weighted_off_center_yy=std::abs(Weighted_CY_Center-0.5);
411  weighted_off_center_yz=std::abs(Weighted_PY_Center-0.5);
412  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)));
413 }
414 
416  int max=std::min(PX_max,PY_max);
417  int min=std::max(PX_min,PY_min);
418  double number_hits_in_overlap_planes=0;
419  for (auto const hit: hits) {
420  if(hit.Plane().val<=max && hit.Plane().val>=min) number_hits_in_overlap_planes++;
421  }
422 
423  return number_hits_in_overlap_planes/hits.size();
424 
425 }
426 
T max(const caf::Proxy< T > &a, T b)
def stdev(lst)
Definition: HandyFuncs.py:286
void StdevCellsPerPlane(HitList const &hits, double &stdev_cells_per_plane_xview, double &stdev_cells_per_plane_yview) const
set< int >::iterator it
virtual bool filter(art::Event &e) override
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
novaddt::Plane const & Plane() const
Definition: DAQHit.h:70
std::vector< DAQHit > HitList
Definition: HitList.h:15
const char * p
Definition: xmltok.h:285
T sqrt(T number)
Definition: d0nt_math.hpp:156
value_type val
Definition: BaseProducts.h:109
DEFINE_ART_MODULE(TestTMapFile)
unsigned distance(const T &t1, const T &t2)
float abs(float number)
Definition: d0nt_math.hpp:39
Identifier for the Y measuring view of the detector (side)
ProductID put(std::unique_ptr< PROD > &&product)
Definition: Event.h:102
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
unsigned long long value_type
Definition: BaseProducts.h:32
void hits()
Definition: readHits.C:15
float Distance(unsigned const &xcell1, unsigned const &ycell1, unsigned const &plane1, unsigned const &xcell2, unsigned const &ycell2, unsigned const &plane2) const
Identifier for the X measuring view of the detector (top)
value_type val
Definition: BaseProducts.h:84
double NumberOfHitsInOverlapPlanesCut(HitList const &hits, float PX_min, float PX_max, float PY_min, float PY_max) const
static float min(const float a, const float b, const float c)
Definition: absgeo.cxx:45
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
Definition: structs.h:12
novaddt::Cell const & Cell() const
Definition: DAQHit.h:71
FastMMTrigger(fhicl::ParameterSet const &p)
bool SurfAssign(DAQHit const &hit) const
T min(const caf::Proxy< T > &a, T b)
Double_t sum
Definition: plot.C:31
Float_t e
Definition: plot.C:35
T max(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:68
int NumberOfSurfaceHits(HitList const &hits) const
enum BeamMode string