SPCluster_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // Class: SPCluster
3 // Module Type: producer
4 // File: SPCluster_module.cc
5 //
6 // Generated at Wed Dec 18 15:41:32 2013 by Zukai Wang using artmod
7 // from cetpkgsupport v1_04_02.
8 ////////////////////////////////////////////////////////////////////////
16 #include "fhiclcpp/ParameterSet.h"
18 
19 #include "Geometry/Geometry.h"
22 
23 #include "NovaDAQConventions/DAQConventions.h"
24 #include "RecoBase/CellHit.h"
25 #include "RecoBase/Cluster.h"
26 
27 #include "TH3.h"
28 #include "TVector3.h"
29 #include "TTree.h"
30 
31 #include <algorithm>
32 #include <utility>
33 #include <cmath>
34 
35 namespace zcl {
36  class SPCluster;
37 }
38 
40 public:
41  explicit SPCluster(fhicl::ParameterSet const & pset);
42  virtual ~SPCluster();
43 
44  void produce(art::Event & evt);
45  void beginJob();
46  void endJob();
47 
48  bool SurfAssign(art::Ptr<rb::CellHit> const & hptr);
49  bool ContaineAssign(art::Ptr<rb::CellHit> const & hptr);
50 
51  //Return the indices of the earliest and latest matching elements in the other view:
52  bool ViewCheck(art::Ptr<rb::CellHit> const &,
53  art::Ptr<rb::CellHit> const &);
54 
55  float Dist(art::Ptr<rb::CellHit> const &,
56  art::Ptr<rb::CellHit> const &);
57 
58  //This function return the index of the hit which earlier than and closest to the time you assigned
59  unsigned IndexLocator(std::vector<art::Ptr<rb::CellHit> > & hit_ptrs,
60  float tns);
61 
62  bool dupCheck(std::vector<art::Ptr<rb::CellHit> > & hit_ptrs,
63  art::Ptr<rb::CellHit>& chit);
64 
65 private:
66  // FHICL Parameters:
68  float _hf_min; //The threshold of hits/distance
69  //One should consider the speed of MM going to the next plane in the other view,
70  //so the timing cut should be looser in accomodating slower monopoles
71  unsigned _sigmaT; //Time resolution (ns) for single hit: same view parameter
72  unsigned _deltaT; //for checking mate hits on the other view
73  unsigned _minHits; //precut on hits number limit
74  float _deltaP_S; //This is used for surface determination and also for view matching determination
75  float _deltaC_S; //For surface determination
76  float _deltaP_C; //This is used for surface determination and also for view matching determination
77  float _deltaC_C; //For surface determination
78  float _sigmaL; // Triangle cut for the max additional ratio allowed
79  float _Xmax;
80  float _Ymax;
81  float _DistMin; //Min distance requirement
82  float _vproMin; //Max beta: 0.01
83  float _vproMax; //Min beta: 4e-5
84  float _df; //Maximum fraction of (plane gap)/(plane range) allowed
85 
86 };
87 
89  fCellHitInput (p.get< std::string >("CellHitInput") ),
90  _hf_min (p.get< float >("hf_min")), //Min hit number density in unit distance(in cell unit)
91  _sigmaT (1500), //100TDC = 1562.5ns
92  _deltaT (4687), //300TDC
93  _minHits (30),
94  _deltaP_S (4),
95  _deltaC_S (10),
96  _deltaP_C (2),
97  _deltaC_C (2),
98  _sigmaL (0.02),
99  _Xmax (p.get<float>("Xmax")),
100  _Ymax (p.get<float>("Ymax")),
101  _DistMin (p.get<float>("DistMin")),
102  _vproMin (p.get<float>("vproMin")),
103  _vproMax (p.get<float>("vproMax")),
104  _df (0.1)
105 {
106  // Call appropriate Produces<>() functions here.
107  produces< std::vector<rb::Cluster> >();
108 }
109 
111 {
112  // Clean up dynamic memory and other resources here.
113 }
114 
116 {
117 }
118 
120 {
123  evt.getByLabel(fCellHitInput, hitcol);
124  // Load the cell hits into a vector for easy use
125  std::vector<art::Ptr<rb::CellHit > > hitlist[2];
126 
127  //This process imitated the ADC filtering process
128  for(unsigned int i = 0; i < hitcol->size(); ++i){
129  art::Ptr<rb::CellHit> hit_ptr(hitcol, i);
130  //Make a selection of hits here: first 2 diblocks
131 
132  if (hit_ptr->ADC() > 50 && hit_ptr->Plane() < _Xmax+1) {
133  unsigned v = 0;
134  if (hit_ptr->View()==geo::kY ) v = 1;
135  hitlist[v].emplace_back(hit_ptr);
136  }
137  }
138 
139 
140  //Sort the hits by time:
141  rb::SortByTime(hitlist[0]);
142  rb::SortByTime(hitlist[1]);
143 
144  //X, Y, T (plane, cell, tns)
145  std::vector<art::Ptr<rb::CellHit> > c_PS[2]; //They will contain all Valid non surface hits on each view
146  std::vector<art::Ptr<rb::CellHit> > s_PS[2]; //They will contain all Valid surface and mate hits on each view
147 
148  std::unique_ptr<std::vector<rb::Cluster> > clustercol(new std::vector<rb::Cluster>);
149 
150  //View Check and doing surface hits collecion
151  for (unsigned v = 0; v!=2; ++v) {
152  for (unsigned i = 0; i!= hitlist[v].size(); ++i ) {
153  art::Ptr<rb::CellHit> const& hit = hitlist[v][i];
154  float t = hit->TNS();
155  unsigned start_id = IndexLocator(hitlist[1-v], t-_deltaT);
156  unsigned stop_id = IndexLocator(hitlist[1-v], t+_deltaT);
157  bool met = false;
158  for (unsigned j = start_id; j<stop_id; ++j ) {
159  art::Ptr<rb::CellHit> const& buddy = hitlist[1-v][j];
160  if (ViewCheck(hit, buddy)) {
161  if (!met && ContaineAssign(hit) ) {
162  c_PS[v].emplace_back(hit);
163  met=true;
164  }
165  else if (SurfAssign(hit) || SurfAssign(buddy)) {
166  s_PS[v].emplace_back(hit);
167  break;
168  }
169  }
170  }
171  }
172  }
173 
174  for (unsigned v=0; v!=2; ++v) {
175  if (c_PS[v].size()<_minHits) {
176  evt.put(std::move(clustercol));
177  return;
178  }
179  if (s_PS[v].size()<2) {
180  evt.put(std::move(clustercol));
181  return;
182  }
183  }
184 
185  for (int v=0; v!=2; ++v) {
186  // std::cout<<"View: "<<v<<std::endl;
187  //If all X view hits have all been reviewed and no cluster has been reconstructed, you can stop here:
188  if (v==1 && clustercol->empty() ) break;
189 
190  //Try to find entry hit on valid surface hits containers
191  for (unsigned i=0; i<s_PS[v].size()-1; ++i) {
192  // bool found = false;
193  art::Ptr<rb::CellHit> hit1 = s_PS[v][i];
194  float ti = hit1->TNS();
195  int pi = hit1->Plane();
196  int ci = hit1->Cell();
197 
198  // std::cout<<"potential Pi: "<<pi<<" Ci: "<<ci<<std::endl;
199 
200 
201  unsigned ini = IndexLocator(c_PS[v], ti-_deltaT);
202  // std::cout<<"ini:"<<ini<<std::endl;
203 
204  if (ini+_minHits+1 > c_PS[v].size() ) break;
205 
206  for (unsigned j=s_PS[v].size()-1; j>i && ini+_minHits+1 < c_PS[v].size(); --j) {
207  art::Ptr<rb::CellHit> hit2 = s_PS[v][j];
208  float tf = hit2->TNS();
209  int pf = hit2->Plane();
210  int cf = hit2->Cell();
211 
212  if (tf +_deltaT < c_PS[v][ini+_minHits]->TNS()) {
213  // std::cout<<"break due to tf too small"<<std::endl;
214  break;
215  }
216  if (tf > c_PS[v][c_PS[v].size()-1]->TNS() +_deltaT ){
217  // std::cout<<"continue due to tf too large"<<std::endl;
218  continue;
219  }
220 
221  unsigned fin = IndexLocator(c_PS[v], tf+_deltaT);
222  float Distance = Dist(hit1, hit2);
223 
224  if (Distance < _DistMin) {
225  continue;
226  }
227  if (fin<ini+Distance*_hf_min) {
228  // std::cout<<"fin:"<<fin<<std::endl;
229  continue;
230  }
231 
232  float vpro = (tf-ti)/Distance;
233  // std::cout<<"vpro: "<<vpro<<std::endl;
234  if (tf > ti+Distance*_vproMax+_sigmaT) {
235  // std::cout<<"tf too large"<<std::endl;
236  continue;
237  }
238  else if (tf+_sigmaT < ti+Distance*_vproMin) {
239  // std::cout<<"tf too small"<<std::endl;
240  continue;
241  }
242 
243  //Here we implemente the CCB(continuity-check and break) mechanisim
244  //First, calculate the time window for checking the previous hit on track:
245  int _dP = abs(pf-pi)*_df;
246  float dt = (tf-ti)*_df;
247 
248  //Fish is coming!
249  unsigned h_counts = 0;
250  float disCut = std::min(Distance*_sigmaL,float(3));
251 
252  unsigned h_cut = Distance*_hf_min;
253  rb::Cluster Cluster;
254 
255  int last_P = pi;
256  float last_T = ti;
257  for (unsigned k=ini; k!=fin; ++k) {
258  art::Ptr<rb::CellHit> Mid = c_PS[v][k];
259  int pk = Mid->Plane();
260  int ck = Mid->Cell();
261  float tk = Mid->TNS();
262 
263  // std::cout<<"k-ini: "<<k-ini<<std::endl;
264  if (tk>last_T+dt+_sigmaT) {
265  // std::cout<<"hits attached: "<<h_counts<<" current Plane: "<<pk<<" curent time: "<<tk-last_T<<" view:"<<v<<std::endl;
266  if (pk> last_P+_dP || pk+_dP<last_P) break;
267  }
268 
269  if (pk > pi && pk > pf) continue;
270  else if (pk < pi && pk < pf) continue;
271  if (ck > ci && ck > cf) continue;
272  else if (ck < ci && ck < cf) continue;
273 
274  float t1 = tk-ti;
275  float t2 = tf-tk;
276  float d1 = Dist(hit1, Mid);
277  float d2 = Dist(hit2, Mid);
278 
279  if (d1+d2 > Distance+disCut) continue;
280 
281  if (d1>d2) {
282  if (t1 < d1*vpro + _sigmaT && t1 > d1*vpro - _sigmaT) {
283  ++h_counts;
284  Cluster.Add(Mid);
285  //refresh the previous hit:
286  last_P= pk;
287  last_T= tk;
288  }
289  }
290  else {
291  if (t2 > d2*vpro - _sigmaT && t2 < d2*vpro + _sigmaT) {
292  ++h_counts;
293  Cluster.Add(Mid);
294  //refresh the previous hit:
295  last_P= pk;
296  last_T= tk;
297  }
298  }
299  }
300 
301  if (h_counts > h_cut ) {
302  //Now delete the hits from the container:c_PS[v]
303  for (unsigned idx =0; idx != Cluster.NCell(); ++idx) {
304  auto it = find(c_PS[v].begin(), c_PS[v].end(), Cluster.Cell(idx));
305  c_PS[v].erase(it);
306  }
307 
308  Cluster.Add(hit1);
309  Cluster.Add(hit2);
310  clustercol->emplace_back(Cluster);
311 
312 
313  }
314  }
315  } // 1 view finished
316  }
317 
318  evt.put(std::move(clustercol));
319 } //end production function
320 
322 
323 }
324 
326 {
327  if ( hptr->Cell() > _Ymax - _deltaC_S ) return true;
328  else if ( hptr->Cell() < _deltaC_S ) return true;
329  else if ( hptr->Plane() < _deltaP_S ) return true;
330  else if ( hptr->Plane() > _Xmax - _deltaP_S ) return true;
331 
332  return false;
333 
334 }
335 
337 {
338  if ( hptr->Cell() > _Ymax - _deltaC_C ) return false;
339  else if ( hptr->Cell() < _deltaC_C ) return false;
340  else if ( hptr->Plane() < _deltaP_C ) return false;
341  else if ( hptr->Plane() > _Xmax - _deltaP_C ) return false;
342 
343  return true;
344 }
345 
346 //Note: Ts < Tm (Looking forwards). Time check is already done in the main loop.
348  art::Ptr<rb::CellHit> const & hptr2)
349 {
350  if (hptr1->TNS() > hptr2->TNS()+_deltaT) return false;
351  if (hptr2->TNS() > hptr1->TNS()+_deltaT) return false;
352  if (hptr1->Plane() > hptr2->Plane() ) {
353  if (hptr1->Plane() - hptr2->Plane() < _deltaP_S ) return true;
354  }
355  else if (hptr2->Plane() - hptr1->Plane() < _deltaP_S ) return true;
356  return false;
357 }
358 
360  art::Ptr<rb::CellHit> const & hptr2)
361 {
362  float dx = ((hptr1->Plane() )- (hptr2->Plane()) )*1.6795;
363  float dy = (hptr1->Cell() ) - (hptr2->Cell());
364  return sqrt(dx*dx+dy*dy);
365 }
366 
367 //It should return the index before T
368 unsigned zcl::SPCluster::IndexLocator(std::vector<art::Ptr<rb::CellHit> > & hit_ptrs,
369  float tns)
370 {
371  unsigned i = 0;
372  unsigned f = hit_ptrs.size()-1;
373 
374  if (tns < (hit_ptrs[0])->TNS()) return i;
375  if (tns > (hit_ptrs[f])->TNS()) return f;
376 
377  while (f>i+1) {
378  unsigned m = (f+i)/2;
379  if (tns >= (hit_ptrs[m])->TNS() ) i = m;
380  else f = m;
381  }
382  return i;
383 }
384 
385 bool zcl::SPCluster::dupCheck(std::vector<art::Ptr<rb::CellHit> > & hit_ptrs,
386  art::Ptr<rb::CellHit> & chit)
387 {
388  if (hit_ptrs.empty()) return false;
389  float t = chit->TNS();
390 
391  unsigned f = hit_ptrs.size()-1;
392  if (t > hit_ptrs[f]->TNS()) return false;
393  unsigned i = 0;
394 
395  while (f>i+5) {
396  unsigned m = (f+i)/2;
397  if (t >= hit_ptrs[m]->TNS() ) i = m;
398  else f = m+1;
399  }
400 
401  for (unsigned idx = i; idx!=f; ++idx) {
402  if (chit==hit_ptrs[idx])
403  return true;
404  }
405 
406  return false;
407 
408 }
409 
SPCluster(fhicl::ParameterSet const &pset)
float TNS() const
Definition: CellHit.h:46
TString fin
Definition: Style.C:24
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
Definition: Cluster.cxx:134
set< int >::iterator it
unsigned short Plane() const
Definition: CellHit.h:39
geo::View_t View() const
Definition: CellHit.h:41
const char * p
Definition: xmltok.h:285
T sqrt(T number)
Definition: d0nt_math.hpp:156
bool ContaineAssign(art::Ptr< rb::CellHit > const &hptr)
A collection of associated CellHits.
Definition: Cluster.h:47
double Distance(double x1, double y1, double x2, double y2)
void abs(TH1 *hist)
void SortByTime(std::vector< art::Ptr< rb::CellHit > > &c)
Sort c in time order (earliest to latest).
Definition: CellHit.cxx:134
DEFINE_ART_MODULE(TestTMapFile)
Timing fit.
Horizontal planes which measure Y.
Definition: PlaneGeo.h:29
virtual void Add(const art::Ptr< rb::CellHit > &cell, double weight=1)
Definition: Cluster.cxx:84
ProductID put(std::unique_ptr< PROD > &&product)
Definition: Event.h:102
unsigned short Cell() const
Definition: CellHit.h:40
Definition: Cand.cxx:23
std::string fCellHitInput
double dy[NP][NC]
double dx[NP][NC]
bool ViewCheck(art::Ptr< rb::CellHit > const &, art::Ptr< rb::CellHit > const &)
int evt
bool SurfAssign(art::Ptr< rb::CellHit > const &hptr)
const double j
Definition: BetheBloch.cxx:29
double t2
art::Ptr< rb::CellHit > Cell(geo::View_t view, unsigned int viewIdx) const
Get the ith cell from view view.
Definition: Cluster.cxx:145
bool dupCheck(std::vector< art::Ptr< rb::CellHit > > &hit_ptrs, art::Ptr< rb::CellHit > &chit)
int16_t ADC(uint32_t i) const
Definition: RawDigit.cxx:58
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
Definition: structs.h:12
void geom(int which=0)
Definition: geom.C:163
void produce(art::Event &evt)
T min(const caf::Proxy< T > &a, T b)
unsigned IndexLocator(std::vector< art::Ptr< rb::CellHit > > &hit_ptrs, float tns)
Encapsulate the geometry of one entire detector (near, far, ndos)
float Dist(art::Ptr< rb::CellHit > const &, art::Ptr< rb::CellHit > const &)
enum BeamMode string