SlowMMTrigger_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // Class: SlowMMTrigger
3 // Module Type: filter
4 // File: SlowMMTrigger_module.cc
5 //
6 // Generated at Mon Aug 26 13:26:30 2013 by Zukai Wang using artmod
7 // from cetpkgsupport v1_02_00.
8 ////////////////////////////////////////////////////////////////////////
9 
10 //ART Framework Include
16 #include "fhiclcpp/ParameterSet.h"
22 
23 // NOvADDT includes
35 
36 #include <algorithm>
37 #include <fstream>
38 #include <time.h>
39 #include <fstream>
40 
41 namespace novaddt {
42  class SlowMMTrigger;
43 }
44 
46 public:
47  explicit SlowMMTrigger(fhicl::ParameterSet const & p);
48  virtual ~SlowMMTrigger();
49 
50  bool filter(art::Event & e) override;
51 
52  //This function should return an int showidng the district of the hit:
53  //0 -> non surface
54  bool SurfAssign(DAQHit const & hptr);
55 
56  //Return the indices of the earliest and latest matching elements in the other view:
57  bool ViewCheck(unsigned p1, TDC::value_type T1,
58  unsigned p2, TDC::value_type T2);
59 
60  //This function return the index of the hit which closest to th time you assigned
61  unsigned IndexLocator(HitList & hitlist,
63 
64 
65  float Dist_(float c1, float c2,
66  float p1, float p2);
67 
68 private:
69  unsigned _prescale;
70 
73 
74  float _hf_min; //The threshold of hits/distance
75 
76  unsigned _deltaT; //Time resolution (TDC) for single hit
77  unsigned _sigmaT;
78  unsigned _minHits; //precut on hits number limit
79  float _deltP; //This is used for surface determination and also for view matching determination
80  float _deltC; //For surface determination
81  float _sigmaL; // Triangle cut for the max additional ratio allowed
82  float _Pmin;
83  float _Pmax;
84  float _Cmax;
85  float _DistMin; //Min distance requirement
86  float _vproMin; //Max beta: 0.01
87  float _vproMax; //Min beta: 4e-5
88  unsigned _trigger_counts;
89  float _df; //CCB mechanism parameter
90  int _CPUticks; //If you consume more than this amount of CPU time, just trigger the lowest safe timewindow
91 };
92 
94  :
95  _prescale (p.get< unsigned >("prescale")),
96  _hitslabel (p.get< std::string >("hits_label")),
97  _instance (p.get< std::string >("instance" )),
98  _hf_min (p.get< float >("hf_min")), //Min hit number density in unit distance
99  _deltaT (128),
100  _sigmaT (64), //Time difference allowence when attaching hits online & index locator
101  _minHits (p.get<unsigned>("minHits")),
102  _deltP (p.get<float>("deltP")),
103  _deltC (p.get<float>("deltC")),
104  _sigmaL (p.get<float>("sigmaL")),
105  _Pmin (p.get<float>("Pmin")),
106  _Pmax (p.get<float>("Pmax")),
107  _Cmax (p.get<float>("Cmax")),
108  _DistMin (p.get<float>("DistMin")),
109  _vproMin (p.get<float>("vproMin")),
110  _vproMax (p.get<float>("vproMax")),
111  _trigger_counts(0),
112  _df (0.1),
113  _CPUticks (p.get<int>("CPUticks"))
114 {
115  // Call appropriate Produces<>() functions here.
116  produces<std::vector<TriggerDecision>>();
117 }
118 
120 {
121  // Clean up dynamic memory and other resources here.
122 }
123 
125 {
126  //Record the starting time:
127  clock_t start_T = clock();
128  bool result = false;
129  std::unique_ptr<std::vector<TriggerDecision>>
130  trigger_decisions(new std::vector<TriggerDecision>);
131 
132  //Get the hits:
134  e.getByLabel(_hitslabel, _instance, hits);
135 
136  HitList hl[2];
137  //Seperate the hits by view, preserve the time order:
138  for (unsigned i = 0; i!= hits->size(); ++i ) {
139  unsigned v = 0;
140  if ( ((*hits)[i]).View().val==daqchannelmap::Y_VIEW ) v = 1;
141  hl[v].emplace_back((*hits)[i]);
142  }
143 
144  // std::cout<<"Xhits: "<<hl[0].size()<<" Yhits: "<<hl[1].size();
145  if (hl[0].size()<_minHits+2 || hl[1].size()<_minHits+2) {
146  e.put(std::move(trigger_decisions));
147  return false;
148  }
149 
150  //Only record X hits
151  HitList c_PS; //They will contain all Valid non surface hits on each view
152  HitList s_PS; //They will contain all Valid surface and mate hits on each view
153 
154  //View check
155  //Goal: doing view check for each X hit, and preserve the time order
156  for (unsigned i = 0; i!= hl[0].size(); ++i ) {
157  DAQHit const& hx = hl[0][i];
158  TDC::value_type Tx = hx.TDC().val;
159  unsigned p1= hx.Plane().val;
160  //Do index locator calculation to locate the indices range:
161  unsigned start_id = IndexLocator(hl[1], Tx-_deltaT);
162  unsigned stop_id = IndexLocator(hl[1], Tx+_deltaT);
163  bool met = false;
164  //Look for your mate, once found, push back and break:
165  for (unsigned j = start_id; j < stop_id; ++j) {
166  DAQHit const& hy = hl[1][j];
167  unsigned p2 = hy.Plane().val;
168  TDC::value_type Ty = hy.TDC().val;
169 
170  if (ViewCheck(p1,Tx,p2,Ty) ) {
171  if (SurfAssign(hx) || SurfAssign(hy)) {
172  s_PS.emplace_back(hx);
173  break;
174  }
175  else if (!met) c_PS.emplace_back(hx);
176  met = true;
177  }
178  }
179  }
180  /*
181  std::cout<<" ViewCheckTicks: "<<clock()-start_T;
182  std::cout<<" SurfXhits: "<<s_PS.size()
183  <<" ContainedXhits: "<<c_PS.size();
184  */
185 
186  unsigned x_con = c_PS.size();
187  if (x_con < _minHits+2 || s_PS.size()<2) {
188  e.put(std::move(trigger_decisions));
189  return false;
190  }
191 
192  //Making possible combinations of entry & exit hits:
193  //Hints: you usually get more y hits than x hits, so look for candidates in X view
194  for (unsigned i=0; i<s_PS.size()-1 && !result; ++i) {
195  DAQHit const& entry = s_PS[i];
196  TDC::value_type Ti = entry.TDC().val;
197  unsigned pi = entry.Plane().val;
198  unsigned ci = entry.Cell().val;
199  unsigned ini = IndexLocator(c_PS, Ti)+1;
200 
201  if (ini+_minHits+1 > x_con ) break; //Don't go to the next hit, you won't get enough hits
202 
203  //Notice: the entry hit and the exit hit may not be in the same view:
204  //Notice: exit hit should be looping from back to prevent double triggering!
205  //According to the chronological order of hits, once j is triggered,
206  //you don't have to go down further, later triggers will be included in the time window!
207  for (unsigned j=s_PS.size()-1; j>i && !result; --j) {
208  DAQHit const& exit = s_PS[j];
209  TDC::value_type Tf = exit.TDC().val;
210 
211  unsigned pf = exit.Plane().val;
212  unsigned cf = exit.Cell().val;
213 
214  if (Tf+_deltaT < c_PS[ini+1+_minHits].TDC().val ) break;
215  else if (Tf > c_PS[x_con-1].TDC().val+_deltaT ) continue;
216 
217  //Set the plane range:
218  unsigned p_max = std::max(pi,pf);
219  unsigned p_min = std::min(pi,pf);
220  //Set the cell range:
221  unsigned c_max = std::max(ci,cf);
222  unsigned c_min = std::min(ci,cf);
223 
224  //If you are not traveling through enough planes, you also get cut:
225  if (p_max < p_min+2*_deltP) continue;
226 
227  float Distance = Dist_(ci,cf,pi,pf);
228  float vpro = (Tf-Ti)/Distance; // Always positive since hits are sorted
229 
230  if (Distance < _DistMin) continue;
231  if (Tf > Ti + unsigned(_vproMax*Distance)+_sigmaT ) continue;
232  if (Tf+_sigmaT < Ti + unsigned(_vproMin*Distance) ) continue;
233 
234  unsigned fin = IndexLocator(c_PS, Tf);
235 
236  if (fin<ini+Distance*_hf_min) continue;
237 
238  //Before looping the mid point one by one, check if you are running out of time,
239  if (clock()-start_T > _CPUticks) {
240  result = true;
241  trigger_decisions->
242  emplace_back(Ti, Tf - Ti,
244 
245  /*
246  std::cout<<" ot: "<<1<<" ticks: "<<clock()-start_T
247  <<" DT: "<<Tf - Ti<<std::endl;
248  */
249  e.put(std::move(trigger_decisions));
250  return result;
251  }
252 
253 
254  //Hmm, here comes the meat of this piece of code:
255  //Loop through the non-boundary hits in the view, in the given range:
256  unsigned h_counts = 0;
257 
258  float disCut = Distance*_sigmaL;
259  unsigned h_cut = Distance*_hf_min;
260 
261  //Here we implement the CCB(continuity-check and break) mechanisim
262  //First, calculate the time window for checking the previous hit on track:
263  unsigned _dP = (p_max-p_min)*_df;
264  if (_dP<2) _dP=2;
265  else if (_dP>4) _dP=4;
266  float dt = (Tf-Ti)*_df;
267  if (dt>256) dt=256; //4us, enough for beta-4 to travel 12 cm!
268  TDC::value_type T_last = Ti;
269  unsigned P_last = pi;
270  unsigned C_last = ci;
271 
272  for (unsigned k=ini; k!=fin;++k) {
273  DAQHit const& midP = c_PS[k];
274  unsigned pk = midP.Plane().val;
275  unsigned ck = midP.Cell().val;
276  TDC::value_type Tk = midP.TDC().val;
277  if (pk>p_max) continue;
278  if (pk<p_min) continue;
279  if (ck>c_max) continue;
280  if (ck<c_min) continue;
281 
282  //Now you can check whether it is going forward:
283  if (pf>pi) {
284  if ( pk> P_last+_dP || pk < P_last){
285  if (Tk-T_last>dt) break;
286  else continue;
287  }
288 
289  if (cf > ci && ck < C_last ) continue;
290  if (cf < ci && ck > C_last ) continue;
291  }
292 
293  else {
294  if (pk+_dP < P_last || pk > P_last ) {
295  if (Tk-T_last>dt) break;
296  else continue;
297  }
298  if (pk+_dP < P_last || pk > P_last || ck > C_last+_deltC || ck+_deltC<C_last) continue;
299  if (cf > ci && ck < C_last ) continue;
300  if (cf < ci && ck > C_last ) continue;
301  }
302 
303  float d1 = Dist_(ci,ck,pi,pk);
304  float d2 = Dist_(cf,ck,pf,pk);
305 
306  if (d1+d2 > Distance+disCut) continue;
307  float t1 = Tk-Ti; // From start to mid
308  float t2 = Tf-Tk; // From mid to end
309 
310  if (d1>d2) {
311  if (t1 < d1*vpro + _sigmaT && t1 + _sigmaT > d1*vpro) {
312  ++h_counts;
313  //refresh the previous hit:
314  P_last = pk;
315  C_last = ck;
316  T_last = Tk;
317  }
318  }
319 
320  else {
321  if (t2 < d2*vpro + _sigmaT && t2 + _sigmaT > d2*vpro ) {
322  ++h_counts;
323  //refresh the previous hit:
324  P_last = pk;
325  C_last = ck;
326  T_last = Tk;
327  }
328  }
329  }
330 
331  if (h_counts>h_cut && h_counts>_minHits-2) {
332  result = true;
333  trigger_decisions->
334  emplace_back(Ti-2*_deltaT, Tf - Ti+4*_deltaT,
336  /*
337  std::cout<<" ot: "<<0<<" ticks: "<<clock()-start_T
338  <<" DT: "<<Tf - Ti+4*_deltaT;
339  */
340  break;
341  }
342  }
343  }
344  /*
345  std::cout<<" ot: "<<0<<" ticks: "<<clock()-start_T
346  <<" DT: "<<-1<<std::endl;
347  */
348  e.put(std::move(trigger_decisions));
349  return result;
350 
351 }
352 
354 {
355  //Check if this is close to surface
356  if ( hit.Cell().val > _Cmax - _deltC ) return true;
357  else if ( hit.Cell().val < _deltC ) return true;
358  else if ( hit.Plane().val < _deltP ) return true;
359  else if ( hit.Plane().val > _Pmax - _deltP ) return true;
360 
361  return false;
362 }
363 
364 //Note: Ts < Tm (Looking forwards). Time check is already done in the main loop.
366  unsigned p2, TDC::value_type T2)
367 {
368  if (T1>T2+_deltaT) return false;
369  if (T2>T1+_deltaT) return false;
370  if (p1 > p2 ) {
371  if (p1-p2 < _deltP) return true;
372  }
373  else if (p2-p1 < _deltP) return true;
374 
375  return false;
376 }
377 
378 //It should return the index before T
381 {
382  unsigned i = 0;
383  unsigned f = hl.size()-1;
384 
385  if (T+_deltaT < hl[i].TDC().val ) return i;
386  if (T> hl[f].TDC().val+_deltaT ) return f;
387 
388  while (f>i+1) {
389  unsigned m = (f+i)/2;
390  if (T >= hl[m].TDC().val ) i = m;
391  else f = m;
392  }
393  return i;
394 }
395 
396 
398  float p1, float p2)
399 {
400  return sqrt((p1-p2)*(p1-p2)*2.82072+(c1-c2)*(c1-c2));
401 }
402 
403 
T max(const caf::Proxy< T > &a, T b)
TString fin
Definition: Style.C:24
value_type val
Definition: BaseProducts.h:34
SlowMMTrigger(fhicl::ParameterSet const &p)
novaddt::Plane const & Plane() const
Definition: DAQHit.h:70
float Dist_(float c1, float c2, float p1, float p2)
novaddt::TDC const & TDC() const
Definition: DAQHit.h:74
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
double Distance(double x1, double y1, double x2, double y2)
DEFINE_ART_MODULE(TestTMapFile)
Identifier for the Y measuring view of the detector (side)
c2
Definition: demo5.py:33
ProductID put(std::unique_ptr< PROD > &&product)
Definition: Event.h:102
unsigned long long value_type
Definition: BaseProducts.h:32
Definition: Cand.cxx:23
void hits()
Definition: readHits.C:15
bool SurfAssign(DAQHit const &hptr)
bool ViewCheck(unsigned p1, TDC::value_type T1, unsigned p2, TDC::value_type T2)
bool filter(art::Event &e) override
const double j
Definition: BetheBloch.cxx:29
value_type val
Definition: BaseProducts.h:84
double t2
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
unsigned IndexLocator(HitList &hitlist, TDC::value_type T)
Definition: structs.h:12
exit(0)
novaddt::Cell const & Cell() const
Definition: DAQHit.h:71
c1
Definition: demo5.py:24
double T
Definition: Xdiff_gwt.C:5
T min(const caf::Proxy< T > &a, T b)
Float_t e
Definition: plot.C:35
def entry(str)
Definition: HTMLTools.py:26
enum BeamMode string