KalmanGeoHelper.cxx
Go to the documentation of this file.
1 //////////////////////////////////////////////////////////////////////////////
2 // \file KalmanGeoHelper.cxx
3 // \brief A class of useful track geometry functions for the KalmanTrack algorithm
4 // \version $Id: KalmanGeoHelper.cxx,v 1.1 2012-11-16 16:57:10 raddatz Exp $
5 // \author Nick Raddatz
6 ////////////////////////////////////////////////////////////////////////////
8 #include "GeometryObjects/Geo.h"
10 #include "NovaDAQConventions/DAQConventions.h"
11 #include "TFile.h"
12 #include <iostream>
13 
14 namespace trk
15 {
16  static std::vector<double> planez;
17  static std::vector<double> planehalfd;
18 
19  static std::vector< std::vector< std::array<double, 3> > > cellxyz;
20 
22  {
23  if(!planez.empty()) return;
24 
25  double posp[3];
26  for(unsigned int i = 0; i < fgeom->NPlanes(); i++){
27  fgeom->Plane(i)->Cell(0)->GetCenter(posp);
28  planez.push_back(posp[2]);
29  planehalfd.push_back(fgeom->Plane(i)->Cell(0)->HalfD());
30  }
31 
32  cellxyz.resize(fgeom->NPlanes());
33 
34  for(unsigned int p = 0; p < fgeom->NPlanes(); p++){
35  const geo::PlaneGeo* geoplane = fgeom->Plane(p);
36  for(unsigned int c = 0; c < geoplane->Ncells(); c++){
37  double xyz[3];
38  fgeom->Plane(p)->Cell(c)->GetCenter(xyz);
39  std::array<double, 3> vxyz;
40  vxyz[0] = xyz[0];
41  vxyz[1] = xyz[1];
42  vxyz[2] = xyz[2];
43  cellxyz[p].push_back(vxyz);
44  }
45  }
46 
47  }
48 
49  //......................................................................
51  {
52  fBadChannels = badChan;
53  SetupHelper();
54  }
55 
56  //......................................................................
58  {
59  }
60 
61  //......................................................................
63  {
65  if(fgeom->DetId() == novadaq::cnv::kFARDET) fd = true;
66  else if(fgeom->DetId() == novadaq::cnv::kNEARDET) nd = true;
67  else if(fgeom->DetId() == novadaq::cnv::kNDOS) ndos = true;
68  else if(fgeom->DetId() == novadaq::cnv::kTESTBEAM) tb = true;
69  else{
70  std::cerr<<"This geometry is not ndos, tb, nd, or fd. This object doesn't know what to do"<<std::endl;
71  }
72  if(fd){
73  fclipperLength = 1.5;
74  fwlow = -700;
75  fwhigh = 700;
76  flengthmax = 6.0;
77  effpar0 = TF1("par0","pol4",fwlow,fwhigh);
78  effpar1 = TF1("par1","pol4",fwlow,fwhigh);
79  effpar2 = TF1("par2","pol5",fwlow,fwhigh);
80  effpar3 = TF1("par3","pol4",fwlow,fwhigh);
81 
82  // low threshold tuning
83  effpar0.SetParameters(0.204074, -2.45217E-6, -3.61656E-7, 4.38289E-11, 4.3173E-13);
84  effpar1.SetParameters(0.126228, 3.77528E-4, 7.85744E-7, -4.65944E-10, -1.0345E-12);
85  effpar2.SetParameters(0.327222, -5.80372E-4, 7.42265E-7, 1.27296E-9, -1.13353E-12, -1.13839E-15);
86  effpar3.SetParameters(1.35419E-2, -8.28103E-5, 8.92098E-8, 5.64673E-11, -6.05079E-14);
87 
88  }
89  else if(nd){
90  fclipperLength = 1.5;
91  fwlow = -200;
92  fwhigh = 200;
93  flengthmax = 6.0;
94  effpar0 = TF1("par0","pol4",fwlow,fwhigh);
95  effpar1 = TF1("par1","pol4",fwlow,fwhigh);
96  effpar2 = TF1("par2","pol5",fwlow,fwhigh);
97  effpar3 = TF1("par3","pol4",fwlow,fwhigh);
98 
99  // low threshold tuning
100  effpar0.SetParameters(0.139861, -5.05093E-5, 2.66785E-7, -3.65295E-10, -1.46979E-11);
101  effpar1.SetParameters(0.525903, 2.00583E-4, -4.07681E-7, 5.28495E-10, 2.01945E-11);
102  effpar2.SetParameters(0.0239177, -6.04344E-5, -1.08969E-7, -3.69116E-10, 2.01116E-11, 5.79678E-14);
103  effpar3.SetParameters(-1.03056E-3, 1.7984E-5, -3.04568E-7, -7.85076E-10, 3.55217E-12);
104 
105  }
106  else if(ndos){
107  fclipperLength = 1.5;
108  }
109  // not sure how the effpars are defined so for test beam just use nd setup
110  else if(tb){
111  fclipperLength = 1.5;
112  fwlow = -130;
113  fwhigh = 130;
114  flengthmax = 6.0;
115  effpar0 = TF1("par0","pol4",fwlow,fwhigh);
116  effpar1 = TF1("par1","pol4",fwlow,fwhigh);
117  effpar2 = TF1("par2","pol5",fwlow,fwhigh);
118  effpar3 = TF1("par3","pol4",fwlow,fwhigh);
119 
120  // low threshold tuning this produces a dist with relatively flat efficiencies as a fn of W
121  effpar0.SetParameters(0.139861, -5.05093E-5, 2.66785E-7, -3.65295E-10, -1.46979E-11);
122  effpar1.SetParameters(0.525903, 2.00583E-4, -4.07681E-7, 5.28495E-10, 2.01945E-11);
123  effpar2.SetParameters(0.0239177, -6.04344E-5, -1.08969E-7, -3.69116E-10, 2.01116E-11, 5.79678E-14);
124  effpar3.SetParameters(-1.03056E-3, 1.7984E-5, -3.04568E-7, -7.85076E-10, 3.55217E-12);
125  }
126 
127  }
128 
129  //......................................................................
130  double KalmanGeoHelper::BestTrackPoint(double z, double v, double halfz, double halfv, std::array<double, 2> & bestPos,
131  double z0,double v0, double dz, double dv)
132  {
133  // bestPos is a vector holding the best possible track hit position in the cell: bestPos[0] = z bestPos[1] = x/y.
134  // output vector the minimum distance between the projected track and the best track hit position in the cell.
135  // This function will take the initial state and z position and try to find the most compatible track hit point with the given hit
136 
137  double slope = dv/dz;
138  // Get the projected track at the cell boundaries.
139  double zlo = z - halfz;
140  double zhi = z + halfz;
141  double vlo = v - halfv;
142  double vhi = v + halfv;
143  double v1 = slope*(zlo-z0)+v0;
144  double v2 = slope*(zhi-z0)+v0;
145  double z1 = (vlo-v0)/slope+z0;
146  double z2 = (vhi-v0)/slope+z0;
147  double deltaz = 0;
148  double deltav = 0;
149  // Make a decision as to how to look for best track point as to avoid dividing by zero
150  if(dz == 0){
151  // either the z is in the z range or else it never will be, either case best position is z0, v
152  bestPos[0] = z0;
153  bestPos[1] = v;
154  // Increment deltav accordingly if actually in the cell
155  if(z >= zlo && z <= zhi){ deltav = vhi - vlo; }
156  } // end of if dz == 0
157  else{
158  // Go through the various cases to see if the track projection intersects the cell hit or not
159  if(slope>0){
160  // check if track enters cell from left side
161  if(v1 >= vlo && v1 <= vhi){
162  // check if track exits on the right side of the cell
163  if(v2 <= vhi){
164  bestPos[0] = z;
165  bestPos[1] = 0.5*v1+0.5*v2;
166  deltaz = zhi-zlo;
167  deltav = v2-v1;
168  }
169  else{
170  // track exits through the top of the cell
171  bestPos[0] = 0.5*zlo+0.5*z2;
172  bestPos[1] = 0.5*v1 + 0.5*vhi;
173  deltaz = z2-zlo;
174  deltav = vhi-v1;
175  }
176  }
177  // check if track enters cell from bottom side
178  else if(v1 < vlo && v2 >= vlo){
179  // check if track exits on right side of cell
180  if(v2 <= vhi){
181  bestPos[0] = 0.5*z1+0.5*zhi;
182  bestPos[1] = 0.5*vlo+0.5*v2;
183  deltaz = zhi-z1;
184  deltav = v2-vlo;
185  }
186  else{
187  // track exits through top of the cell
188  bestPos[0] = 0.5*z1 + 0.5*z2;
189  bestPos[1] = v;
190  deltaz = z2-z1;
191  deltav = vhi-vlo;
192  }
193  }
194  else{
195  // track does not go through cell
196  // check if track is above or below cell
197  if(v1 < vlo){
198  // closest point is bottom right corner of the cell
199  bestPos[0] = zhi;
200  bestPos[1] = vlo;
201  }
202  else{
203  // closest point is top left corner of the cell
204  bestPos[0] = zlo;
205  bestPos[1] = vhi;
206  }
207  }
208  }
209  else{
210  // check if track enters cell from left side
211  if(v1 >= vlo && v1<= vhi){
212  // check if the track exits on the right side of the cell
213  if(v2 >= vlo){
214  bestPos[0] = z;
215  bestPos[1] = 0.5*v1+0.5*v2;
216  deltaz = zhi-zlo;
217  deltav = v2-v1;
218  }
219  else{
220  // track exits through the bottom of the cell
221  bestPos[0] = 0.5*zlo+0.5*z1;
222  bestPos[1] = 0.5*v1+0.5*vlo;
223  deltaz = zlo-z1;
224  deltav = v1-vlo;
225  }
226  }
227  // check if the track enters the cell from the top side
228  else if(v1 >= vhi && v2 <= vhi){
229  // check if the track exits on the right side of the cell
230  if(v2 >= vlo){
231  bestPos[0] = 0.5*z2+0.5*zhi;
232  bestPos[1] = 0.5*vhi+0.5*v2;
233  deltaz = z2-zhi;
234  deltav = vhi-v2;
235  }
236  else{
237  // track exits through the bottom of the cell
238  bestPos[0] = 0.5*z2+0.5*z1;
239  bestPos[1] = v;
240  deltaz = z2-z1;
241  deltav = vhi-vlo;
242  }
243  }
244  else{
245  // track does not go through the cell
246  // check if the track is above or below the cell
247  if(v1 > vhi){
248  // closest point is top right corner of the cell
249  bestPos[0] = zhi;
250  bestPos[1] = vhi;
251  }
252  else{
253  // closest point is the bottom left corner of the cell
254  bestPos[0] = zlo;
255  bestPos[1] = vlo;
256  }
257  // if slope is truly zero than pick midpoint of cell instead of corner points for best z position
258  if(slope == 0){ bestPos[0] = z; }
259  }
260  }
261  } // end of else for if dz == 0
262  return (deltaz*deltaz+deltav*deltav);
263  }
264 
265  //......................................................................
266 
267  int KalmanGeoHelper::CountMissedCellsOnLine(geo::View_t view,double vi, double zi, double vf, double zf, int maxCells)
268  {
269  static double fDetHalfWidth = fgeom->DetHalfWidth();
270  static double fDetHalfHeight = fgeom->DetHalfHeight();
271  static double fDetLength = fgeom->DetLength();
272  static int fplaneMax = fgeom->NPlanes();
273 
274  fillprivategeo();
275 
276  // if we are in the nd or ndos nead to worry about muon catcher which has different spacing between planes and distance
277  if(nd || ndos){
278  // take the ratio of number of planes to distance from the first 100 planes otherwise will get screwed up by the muon catcher
279  fplaneMax = 100;
280  fDetLength = planez[100];
281  }
282  // Otherwise we would count when the ray was outside the detector in the other view.
283  TVector3 r1(vi,0,zi);
284  TVector3 r2(vf,0,zf);
285  if(view == geo::kY){
286  r1.SetXYZ(0,vi,zi);
287  r2.SetXYZ(0,vf,zf);
288  }
289  geo::ClampRayToDetectorVolume(&r1, &r2, &(*fgeom));
290  int nummiss = 0;
291  if(r1.Z() > r2.Z()) std::swap(r1, r2);
292 
293  double z1 = r1.Z();
294  double z2 = r2.Z();
295 
296  const unsigned int planeMax = fgeom->NPlanes();
297 
298  // find the approximate start plane
299  unsigned int iplane = 0;
300  int pl = (int) z1/(fDetLength/((double)fplaneMax));
301  if(pl < 0){iplane = 0; }
302  else if (pl >= (int)planeMax){iplane = planeMax-1;}
303  else iplane = pl;
304  // check that our initial plane guess is correct
305  // ask if the current position is in this plane
306  if(z1<(planez[iplane]-planehalfd[iplane])){
307  bool foundInitPlane = false;
308  // Did not get the right plane location. Check below if position is lower than the estimated plane center position
309  // only need to worry about geusses that are too high in plane number otherwise accept the estimated guess as the start plane
310  int currTestPlane = iplane-1;
311  while(currTestPlane >= 0 && !foundInitPlane){
312  if(z1>(planez[currTestPlane]-planehalfd[currTestPlane])){
313  iplane = currTestPlane;
314  foundInitPlane = true;
315  }
316  else{
317  if(currTestPlane == 0){iplane = 0;}
318  --currTestPlane;
319  }
320  }
321  }
322 
323  const double v1 = r1[view];
324  const double v2 = r2[view];
325 
326  const double denom = z2-z1;
327  const double dvdz = denom ? (v2-v1)/denom : 0;
328 
329  if(denom == 0 && r2[view] < r1[view]){
330  std::swap(r1, r2);
331  z1 = r1.Z();
332  z2 = r2.Z();
333  }
334 
335  for(unsigned int plane = iplane; plane < planeMax; ++plane){
336  const geo::PlaneGeo* geoplane = fgeom->Plane(plane);
337 
338  // Only look at the planes in the specified view
339  if(geoplane->View() != view) continue;
340  // Assuming here perfect alignment and identical cells across plane
341  const double zc = planez[plane];
342  const double dz = planehalfd[plane];
343 
344  if(z1 > zc+dz) continue;
345  // Plane positions increase monotonically
346  if(z2 < zc-dz) break;
347 
348  // Position entering the cell or starting
349  const double zin = std::max(zc-dz, z1);
350  double vlo = v1+(zin-z1)*dvdz;
351  // Position exiting or ending
352  const double zout = std::min(zc+dz, z2);
353  double vhi = v1+(zout-z1)*dvdz;
354 
355  // if slope is infinite set vlo and vhi as appropriately
356  if(denom == 0){
357  vlo = v1;
358  vhi = v2;
359  }
360 
361  // Ensure they're from low to high instead
362  if(vlo > vhi) std::swap(vlo, vhi);
363 
364  //create an approximate map between position and cell number for cells in the xz view.
365  unsigned int icell = 0;
366  int cell = (int) (vlo+fDetHalfWidth)/(2*fDetHalfWidth/((double) geoplane->Ncells()));
367  int viewind = 0;
368  if(view == geo::kY){
369  cell = (int) (vlo+fDetHalfHeight)/(2*fDetHalfHeight/((double) geoplane->Ncells()));
370  viewind = 1;
371  }
372  if(cell < 0){icell = 0; }
373  else if (cell >= (int) geoplane->Ncells()){icell = geoplane->Ncells()-1;}
374  else icell = cell;
375  // check that our initial cell guess is correct
376  const std::array<double, 3> & pos = cellxyz[plane][icell];
377  double dp = geoplane->Cell(icell)->HalfW();
378  if(vlo<(pos[viewind]-dp)){
379  bool foundInitCell = false;
380  // Did not get the right cell location. Check below if position is lower than the estimated cell center position
381  int currTestCell = icell-1;
382  while(currTestCell >= 0 && !foundInitCell){
383  const std::array<double, 3> & ps = cellxyz[plane][currTestCell];
384  dp = geoplane->Cell(currTestCell)->HalfW();
385  // ask if the current position is in this cell
386  if( vlo>(ps[viewind]-dp)){
387  icell = currTestCell;
388  foundInitCell = true;
389  }
390  else{
391  if(currTestCell == 0){icell = 0;}
392  --currTestCell;
393  }
394  }
395  }
396 
397  for(unsigned int jcell = icell; jcell < geoplane->Ncells();++jcell){
398  const std::array<double, 3> & xyz = cellxyz[plane][jcell];
399  const double vc = xyz[view];
400  const double dv = geoplane->Cell(jcell)->HalfW();
401 
402  if(vlo > vc+dv) continue;
403  // Cell positions increase monotonically
404  if(vhi < vc-dv) break;
405  if(fBadChannels){
406  // check if this is a bad channel
407  if(!fbadc->IsBad(plane,jcell) && !fbadc->IsLowEfficiency(plane,jcell)){ ++nummiss; }
408  }
409  else{ ++nummiss; }
410  if(nummiss > maxCells) return nummiss;
411 
412  } // end for jcell
413  } // end for plane
414  return nummiss;
415  } // end of CountMissedCellsOnLine function
416 
417 
418  //......................................................................
419  int KalmanGeoHelper::CountMissedCellsOnLine(geo::View_t view,double vi, double zi, double vf, double zf,std::vector<geo::OfflineChan> &hitsonline, std::vector<double> &distance)
420  {
421  static double fDetHalfWidth = fgeom->DetHalfWidth();
422  static double fDetHalfHeight = fgeom->DetHalfHeight();
423  static double fDetLength = fgeom->DetLength();
424  static int fplaneMax = fgeom->NPlanes();
425 
426  fillprivategeo();
427 
428  // if we are in the nd or ndos nead to worry about muon catcher which has different spacing between planes and distance
429  if(nd || ndos){
430  // take the ratio of number of planes to distance from the first 100 planes otherwise muon catcher will screw this up
431  fplaneMax = 100;
432  fDetLength = planez[100];
433  }
434  // Otherwise we would count when the ray was outside the detector in the other view.
435  TVector3 r1(vi,0,zi);
436  TVector3 r2(vf,0,zf);
437  if(view == geo::kY){
438  r1.SetXYZ(0,vi,zi);
439  r2.SetXYZ(0,vf,zf);
440  }
441  geo::ClampRayToDetectorVolume(&r1, &r2, &(*fgeom));
442  int nummiss = 0;
443  if(r1.Z() > r2.Z()) std::swap(r1, r2);
444 
445  double z1 = r1.Z();
446  double z2 = r2.Z();
447 
448  static const unsigned int planeMax = fgeom->NPlanes();
449 
450  // find the approximate start plane
451  unsigned int iplane = 0;
452  int pl = (int) z1/(fDetLength/((double)fplaneMax));
453  if(pl < 0){iplane = 0; }
454  else if (pl >= (int)planeMax){iplane = planeMax-1;}
455  else iplane = pl;
456  // check that our initial plane guess is correct
457  // ask if the current position is in this plane
458  if(z1<(planez[iplane]-planehalfd[iplane])){
459  bool foundInitPlane = false;
460  // Did not get the right plane location. Check below if position is lower than the estimated plane center position
461  // only need to worry about geusses that are too high in plane number otherwise accept the estimated guess as the start plane
462  int currTestPlane = iplane-1;
463  while(currTestPlane >= 0 && !foundInitPlane){
464  if(z1>(planez[currTestPlane]-planehalfd[currTestPlane])){
465  iplane = currTestPlane;
466  foundInitPlane = true;
467  }
468  else{
469  if(currTestPlane == 0){iplane = 0;}
470  --currTestPlane;
471  }
472  }
473  }
474 
475  const double v1 = r1[view];
476  const double v2 = r2[view];
477 
478  const double denom = z2-z1;
479  const double dvdz = denom ? (v2-v1)/denom : 0;
480 
481  if(denom == 0 && r2[view] < r1[view]){
482  std::swap(r1, r2);
483  z1 = r1.Z();
484  z2 = r2.Z();
485  }
486 
487  for(unsigned int plane = iplane; plane < planeMax; ++plane){
488  const geo::PlaneGeo* geoplane = fgeom->Plane(plane);
489 
490  // Only look at the planes in the specified view
491  if(geoplane->View() != view) continue;
492  // Assuming here perfect alignment and identical cells across plane
493  const double zc = planez[plane];
494  const double dz = planehalfd[plane];
495 
496  if(z1 > zc+dz) continue;
497  // Plane positions increase monotonically
498  if(z2 < zc-dz) break;
499 
500  // Position entering the cell or starting
501  const double zin = std::max(zc-dz, z1);
502  double vlo = v1+(zin-z1)*dvdz;
503  // Position exiting or ending
504  const double zout = std::min(zc+dz, z2);
505  double vhi = v1+(zout-z1)*dvdz;
506 
507  // if slope is infinite set vlo and vhi as appropriately
508  if(denom == 0){
509  vlo = v1;
510  vhi = v2;
511  }
512 
513  // Ensure they're from low to high instead
514  if(vlo > vhi) std::swap(vlo, vhi);
515 
516  //create an approximate map between position and cell number for cells in the xz view.
517  unsigned int icell = 0;
518  int cell = (int) (vlo+fDetHalfWidth)/(2*fDetHalfWidth/((double) geoplane->Ncells()));
519  int viewind = 0;
520  if(view == geo::kY){
521  cell = (int) (vlo+fDetHalfHeight)/(2*fDetHalfHeight/((double) geoplane->Ncells()));
522  viewind = 1;
523  }
524  if(cell < 0){icell = 0; }
525  else if(cell >= (int) geoplane->Ncells()){icell = geoplane->Ncells()-1;}
526  else{ icell = cell; }
527  // check that our initial cell guess is correct
528  const std::array<double, 3> & pos = cellxyz[plane][icell];
529  double dp = geoplane->Cell(icell)->HalfW();
530  if(vlo<(pos[viewind]-dp)){
531  bool foundInitCell = false;
532  // Did not get the right cell location. Check below if position is lower than the estimated cell center position
533  int currTestCell = icell-1;
534  while(currTestCell >= 0 && !foundInitCell){
535  const std::array<double, 3> & ps = cellxyz[plane][currTestCell];
536  dp = geoplane->Cell(currTestCell)->HalfW();
537  // ask if the current position is in this cell
538  if( vlo>(ps[viewind]-dp)){
539  icell = currTestCell;
540  foundInitCell = true;
541  }
542  else{
543  if(currTestCell == 0){icell = 0;}
544  --currTestCell;
545  }
546  }
547  }
548 
549  for(unsigned int jcell = icell; jcell < geoplane->Ncells();++jcell){
550  const std::array<double, 3> & xyz = cellxyz[plane][jcell];
551  const double vc = xyz[view];
552  const double dv = geoplane->Cell(jcell)->HalfW();
553 
554  if(vlo > vc+dv) continue;
555  // Cell positions increase monotonically
556  if(vhi < vc-dv) break;
557  // get the distance in the cell
558  std::array<double, 2> bestPos;
559  double dist2 = BestTrackPoint(xyz[2],vc,dz,dv,bestPos,z1,v1,z2-z1,v2-v1);
560  ++nummiss;
561  hitsonline.push_back(geo::OfflineChan(plane,jcell));
562  distance.push_back(sqrt(dist2));
563  } // end for jcell
564  } // end for plane
565  return nummiss;
566  } // end of CountMissedCellsOnLine function
567 
568  //......................................................................
570  {
571  // if distance in cell is less clipper length, clipper hit return 1.0, unreliable hit efficiency
572  if(dist < fclipperLength){ return 1.0; }
573 
574  // define the length in cell to use
575  double length = dist;
576  if(length > flengthmax){ length = flengthmax; }
577  // define what cell depth to use
578  double depth = w;
579  if(depth < fwlow){ depth = fwlow; }
580  if(depth > fwhigh){ depth = fwhigh; }
581 
582  double effhit = Efficiency(depth,length);
583 
584  if(effhit > 1){ effhit = 1.0; }
585  else if(effhit<0){ effhit = 0.0; }
586  return 1.0-effhit;
587 
588  } // end of CalcMissFrac
589 
590 
591  //......................................................................
592  double KalmanGeoHelper::Efficiency(double w, double length)
593  {
594  double eff = 0.0;
595 
596  if(fd){
597  double llow = 3.4;
598  double lhigh = 3.75;
599 
600  if( length < llow){
601  eff = effpar0.Eval(w)*length + effpar1.Eval(w);
602  }
603  else if( length < lhigh){
604  eff = effpar2.Eval(w)*length + (effpar0.Eval(w)-effpar2.Eval(w))*llow + effpar1.Eval(w);
605  }
606  else{
607  eff = effpar3.Eval(w)*length + (effpar2.Eval(w)-effpar3.Eval(w))*lhigh +
608  (effpar0.Eval(w)-effpar2.Eval(w))*llow + effpar1.Eval(w);
609  }
610  }
611  else if(nd){
612  double llow = 3.0;
613  double lhigh = 3.75;
614 
615  if( length < llow){
616  eff = effpar0.Eval(w)*length + effpar1.Eval(w);
617  }
618  else if( length < lhigh){
619  eff = effpar2.Eval(w)*length + (effpar0.Eval(w)-effpar2.Eval(w))*llow + effpar1.Eval(w);
620  }
621  else{
622  eff = effpar3.Eval(w)*length + (effpar2.Eval(w)-effpar3.Eval(w))*lhigh +
623  (effpar0.Eval(w)-effpar2.Eval(w))*llow + effpar1.Eval(w);
624  }
625  }
626  else if(tb){
627  double llow = 3.0;
628  double lhigh = 3.75;
629 
630  if( length < llow){
631  eff = effpar0.Eval(w)*length + effpar1.Eval(w);
632  }
633  else if( length < lhigh){
634  eff = effpar2.Eval(w)*length + (effpar0.Eval(w)-effpar2.Eval(w))*llow + effpar1.Eval(w);
635  }
636  else{
637  eff = effpar3.Eval(w)*length + (effpar2.Eval(w)-effpar3.Eval(w))*lhigh +
638  (effpar0.Eval(w)-effpar2.Eval(w))*llow + effpar1.Eval(w);
639  }
640  }
641 
642  else{ eff = 1.0; }
643  // unless clipper, ndos shouldn't miss hits (need to check for bad channel definition of efficiency/badchannels outside of this function)
644 
645  return eff;
646  }
647 
648 
649  //......................................................................
650  int KalmanGeoHelper::MatchTrajectoryToPlane(TVector3 traj, int minplane, int maxplane)
651  {
652  int plane = minplane;
653  int mintrialplane = minplane;
654  int maxtrialplane = maxplane;
655  if(minplane < 0){ mintrialplane = 0; }
656  if(maxplane > fnmaxplane-1){ maxtrialplane = fnmaxplane-1; }
657 
658  // Try finding what plane/cell this point is in
659  const geo::CellUniqueId cid = fgeom->CellId(traj.X(),traj.Y(),traj.Z(),1,1,0.1,2);
660 
661  // if found the cell unique id use it; otherwise backout plane from individual plane width checks
662  if(cid){ fgeom->IdToPlane(cid, &plane); }
663  else{
664  for(int ip = maxtrialplane; ip >= mintrialplane; --ip){
665  int trialPlane = ip;
666 
667  const geo::PlaneGeo* geoplane = fgeom->Plane(trialPlane);
668  // get the estimated cell in this plane
669  int mintrialcell = 0;
670  int maxtrialcell = geoplane->Ncells();
671  double v = traj.X();
672  int trialCell = (int) (v+fgeom->DetHalfWidth())/(2*fgeom->DetHalfWidth()/((double) (geoplane->Ncells()-1)));
673  int viewind = 0;
674  if(geoplane->View() == geo::kY){
675  v = traj.Y();
676  trialCell = (int) (v+fgeom->DetHalfHeight())/(2*fgeom->DetHalfHeight()/((double) (geoplane->Ncells()-1)));
677  viewind = 1;
678  }
679 
680  if(trialCell < 0){trialCell = 0; }
681  else if (trialCell >= (int) geoplane->Ncells()){trialCell = geoplane->Ncells()-1;}
682 
683  // check that our initial cell guess is correct
684  const std::array<double, 3> & trialCenter = cellxyz[ip][trialCell];
685  double dp = geoplane->Cell(trialCell)->HalfW();
686 
687  if(v < (trialCenter[viewind]-dp)){
688  // trial cell is too high, try iterating down to cell 0 until find right cell
689  maxtrialcell = trialCell;
690  for(int ic = maxtrialcell; ic >= mintrialcell; --ic){
691  const std::array<double, 3> & tc = cellxyz[trialPlane][ic];
692  dp = geoplane->Cell(ic)->HalfW();
693  trialCell = ic;
694  if(v > (tc[viewind]-dp)){ break; }
695  }
696  }
697  else if(v > (trialCenter[viewind]+dp)){
698  // trial cell is too low, try iterating up to Ncells in plane until find right cell
699  mintrialcell = trialCell;
700  for(int ic = mintrialcell; ic <= maxtrialcell; ++ic){
701  const std::array<double, 3> & tc = cellxyz[trialPlane][ic];
702  dp = geoplane->Cell(ic)->HalfW();
703  trialCell = ic;
704  if(v > (tc[viewind]-dp)){ break; }
705  }
706  }
707 
708  const std::array<double, 3> & tc = cellxyz[trialPlane][trialCell];
709  double trialPlaneWidth = fgeom->Plane(trialPlane)->Cell(trialCell)->HalfD();
710  // starting from max plane first so if the z is larger than the smallest z
711  // on the plane then this is the correct plane
712  if(traj.Z() > (tc[2]-trialPlaneWidth)){
713  plane = trialPlane;
714  break;
715  }
716  }
717  }
718 
719  return plane;
720  } // end of MatchTrajectoryToPlane
721 
722  //......................................................................
723  int KalmanGeoHelper::MatchTrajectoryToPlaneInView(TVector3 traj, geo::View_t view, int minplane, int maxplane)
724  {
725 
726  assert(fgeom->Plane(minplane)->View() == view);
727  assert(fgeom->Plane(maxplane)->View() == view);
728  int plane = MatchTrajectoryToPlane(traj,minplane,maxplane);
729  geo::View_t trialview = fgeom->Plane(plane)->View();
730 
731  fillprivategeo();
732 
733  if(view == geo::kXorY || trialview == view){ return plane; }
734  else{
735  int lowtrialplane = fgeom->NextPlaneOtherView(plane,-1);
736  int hightrialplane = fgeom->NextPlaneOtherView(plane,+1);
737  if(lowtrialplane == geo::kPLANE_NOT_FOUND){
738  if(hightrialplane == geo::kPLANE_NOT_FOUND){ return minplane; }
739  else{ return hightrialplane; }
740  }
741  else if(hightrialplane == geo::kPLANE_NOT_FOUND){
742  if(lowtrialplane >= minplane){ return lowtrialplane; }
743  else{ return minplane; }
744  }
745  else{
746  // check whether the high or low plane is the closest in z
747  if(fabs(planez[lowtrialplane]+planehalfd[lowtrialplane] - traj.Z()) <=
748  fabs(planez[hightrialplane]-planehalfd[hightrialplane] - traj.Z())){ return lowtrialplane; }
749  else{ return hightrialplane; }
750  }
751  }
752  return minplane;
753 
754  } // end of MatchTrajectoryToPlaneInView
755 
756 } // end of trk namespace
This is the version 2, works for multi-track event.
T max(const caf::Proxy< T > &a, T b)
int CountMissedCellsOnLine(geo::View_t view, double vi, double zi, double vf, double zf, int maxCells)
Variation of geo::CountCellsOnlineFast function specific to KalmanTrack algorithm.
KalmanGeoHelper(bool badChan)
double HalfD() const
Definition: CellGeo.cxx:205
void GetCenter(double *xyz, double localz=0.0) const
Definition: CellGeo.cxx:159
static std::vector< double > planez
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
double Efficiency(double w, double length)
X or Y views.
Definition: PlaneGeo.h:30
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
const PlaneGeo * IdToPlane(const CellUniqueId &id, int *iplane) const
const CellGeo * Cell(int icell) const
Definition: PlaneGeo.h:48
const char * p
Definition: xmltok.h:285
T sqrt(T number)
Definition: d0nt_math.hpp:156
static std::vector< double > planehalfd
unsigned int Ncells() const
Number of cells in this plane.
Definition: PlaneGeo.h:43
double HalfW() const
Definition: CellGeo.cxx:191
double DetLength() const
OStream cerr
Definition: OStream.cxx:7
TString ip
Definition: loadincs.C:5
bool ClampRayToDetectorVolume(TVector3 *p0, TVector3 *p1, const GeometryBase *geom)
If either endpoint is outside the detector move it to the edge.
Definition: Geo.cxx:640
const PlaneGeo * Plane(unsigned int i) const
unsigned distance(const T &t1, const T &t2)
Horizontal planes which measure Y.
Definition: PlaneGeo.h:29
double dist
Definition: runWimpSim.h:113
void swap(art::HLTGlobalStatus &lhs, art::HLTGlobalStatus &rhs)
Track finder for cosmic rays.
View_t View() const
Which coordinate does this plane measure.
Definition: PlaneGeo.h:53
Far Detector at Ash River, MN.
length
Definition: demo0.py:21
double CalcMissFrac(double dist, double w, geo::View_t view)
Estimation of the mip missed hit fraction as a function of the dist through a cell and cell depth w...
double dz[NP][NC]
Float_t E
Definition: plot.C:20
Prototype Near Detector on the surface at FNAL.
int MatchTrajectoryToPlaneInView(TVector3 traj, geo::View_t view, int minplane, int maxplane)
Variation of MatchTrajectoryToPlane that limits the matching plane to be in the input view...
novadaq::cnv::DetId DetId() const
Prefer ds::DetectorService::DetId() instead.
Definition: GeometryBase.h:243
Geometry information for a single readout plane.
Definition: PlaneGeo.h:36
Near Detector in the NuMI cavern.
const CellUniqueId CellId(const double &x, const double &y, const double &z, double dxds=0., double dyds=0., double dzds=1., double step=0.01) const
Collect Geo headers and supply basic geometry functions.
double DetHalfHeight() const
art::ServiceHandle< chaninfo::BadChanList > fbadc
z
Definition: test.py:28
unsigned long long int CellUniqueId
Definition: CellUniqueId.h:15
double DetHalfWidth() const
A (plane, cell) pair.
Definition: OfflineChan.h:17
art::ServiceHandle< geo::Geometry > fgeom
static std::vector< std::vector< std::array< double, 3 > > > cellxyz
assert(nhit_max >=nhit_nbins)
double BestTrackPoint(double z, double v, double halfz, double halfv, std::array< double, 2 > &bestPos, double z0, double v0, double dz, double dv)
Estimate the position on a line within a box.
unsigned int NPlanes() const
T min(const caf::Proxy< T > &a, T b)
bool IsBad(int plane, int cell)
Float_t w
Definition: plot.C:20
const unsigned int NextPlaneOtherView(unsigned int p, int d=+1) const
bool IsLowEfficiency(int plane, int cell)
int MatchTrajectoryToPlane(TVector3 traj, int minplane, int maxplane)
Finds the plane best matched to a trajectory point This method takes a trajectory point and finds the...